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Abstract. A central issue in the theory of astrophysical compact objects and heavy 
ion reactions at intermediate and relativistic energies is the Nuclear Equation of State 
(EoS). On one hand, the large and expanding set of experimental and observational 
data is expected to constrain the behaviour of the nuclear EoS, especially at density 
above saturation, where it is directly linked to fundamental processes which can occur 
in dense matter. On the other hand, theoretical predictions for the EoS at high density 
can be challenged by the phenomenological findings. In this topical review paper we 
present the many-body theory of nuclear matter as developed along different years and 
with different methods. Only nucleonic degrees of freedom are considered. We compare 
the different methods at formal level, as well as the final EoS calculated within each 
one of the considered many-body schemes. The outcome of this analysis should help 
in restricting the uncertainty of the theoretical predictions for the nuclear EoS. 
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1. Introduction 

The knowledge of the nuclear Equation of State (EoS) is one of the fundamental goals in 
nuclear physics which has not yet been achieved. The possibility to extract information 
on the nuclear EoS, in particular at high baryon density is restricted to two fields 
of research. The interplay between the theory and the observations of astrophysical 
compact objects is of great relevance in constraining the nuclear EoS. The enormous 
work that has been developing since the last two decades on the study of heavy ion 
reactions at intermediate and relativistic energies is the other pillar on which one can 
hope two build a reasonable model of the nuclear EoS. On the other hand, theoretical 
predictions of the EoS are essential for modeling heavy ion collisions, at intermediate 
and relativistic energies, and the structure of neutron stars, supernova explosions, binary 
collisions of compact stellar objects and their interactions with black holes. In the 
astrophysical context the dynamics is slow enough and the size scale large enough to 
ensure the local equilibrium of nuclear matter, i.e. hydrodynamics can be applied, and 
therefore the very concept of EoS is extremely useful. On the contrary, in nuclear 
collisions the time scale is the typical one for nuclear processes and the size of the 
system is only one order of magnitude larger than the interaction range or possibly of 
the particle mean free path. The physical conditions in the two contexts are therefore 
quite different. Despite that, by a careful analysis of experimental data on heavy ion 
collisions and astrophysical observations it is possible to connect the two realms of 
phenomena which involve nuclear processes at fundamental level, and the EoS provides 
the crucial concept to establish this link. 

From the theoretical point of view the microscopic theory of nuclear matter has 
a long history and impressive progress has been made along the years. In this topical 
review paper we will first review the many-body theory of nuclear matter and compare 
the predictions of different approaches, Sec. 2-8. In Sections 9-10 possible hints from 
astrophysical observations and heavy ion reactions on the nuclear EoS will be critically 
reviewed, with emphasis on the connections that can be established between the two 
fields. 

2. Many-body theory of the EoS. 

The many-body theory of nuclear matter, where only nucleonic degrees of freedom are 
considered, has developed since several decades along different lines and methods. We 
summarize the most recent results in this field and compare the different methods at 
formal level, as well as the final EoS calculated within each one of the considered many- 
body schemes. The outcome of this analysis should help in restricting the uncertainty 
of the theoretical predictions for the nuclear EoS. 

Within the non-relativistic approach the main microscopic methods are the Bethe- 
Brueckner-Goldstone (BBG) approach and the variational method (VM). The Bethe- 
Brueckner-Goldstone is a general many-body method particularly suited for nuclear 
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systems. It has been extensively applied to homogeneous nuclear matter since many 
years and it has been presented in several review articles and textbooks. For a 
pedagogical review see Baldo (1999), where a short historical introduction and extended 
references can be found. Here we restrict the presentation to the basic structure of the 
method, but we will go to some detail in order to prepare the material needed for a 
formal comparison with other methods. We follow closely the presentation of Baldo 
(1999), at least for the more elementary parts. 

Let us suppose for the moment that only a two-body interaction is present. Then 
the Hamiltonian can be written 

h 2 k 2 1 

H = H + Hi = ^ — — ata k + - {k x k 2 \v \k 3 k 4 )a k la k la k4 a k3 , (1) 
k Zm 1 {h} 

where the operators a} (a) are the usual creation (annihilation) operators. The state 
label {k} includes both the three-momentum k and the spin-isospin variables a, r of the 
single particle state. As usual we will represent the interaction matrix elements as in 
Fig. 1, where a particle (hole) state is represented by a line with an up (down) arrow. 

To be definite, for a purely central local interaction the matrix elements read, in 
general 

(a/3H 7 <J) = J rf 3 r 1 rf 3 r 2 0;(r 1 )0;(r 2 )t;(r 1 - r 2 )0 7 ( ri )0 5 (r 2 ) . (2) 

where the 's are the single particle wave functions. The graph of Fig. 1 can be 
interpreted as an interaction of the particle k% and the hole A; 4 which scatter after the 
interaction to ki and k 2 respectively. Incoming arrows in Fig. 1 correspond to states 
appearing on the right of the matrix elements, while the first (second) dot indicates 
their position in the two-body state. 

These graphs for the matrix elements of v are the building blocks for the more 
complete graphs representing the energy perturbation expansion. 

The starting point of the perturbation expansion is the Gell-Mann and Low theorem 
(Gell-Mann and Low 1951). The theorem is quite general and applies to all systems 
which possess a non-degenerate ground state (with a finite energy). If we call \ipo) the 
ground state of the full Hamiltonian H, the theorem states that it can be obtained from 
the ground state \<f> ) of the unperturbed Hamiltonian H (in our case the free Fermi 




Figure 1. Graphical representation of the NN interaction matrix element. 
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gas ground state) by a procedure usually called the adiabatic "switching on" of the 
interaction 

which entails the normalization {4>o\ipo) — 1- 111 Eq. (3), U^ £ \t) is the evolution operator, 
in the interaction picture, from the generic time t to the time t — of the modified 
Hamiltonian 

H\t) = H + e-^H, , (4) 

where e > 0. Equation (4) implies that the Hamiltonian H e coincides with H in the 
limit t — > — oo and with H at t — and that the interaction is switched on following 
an infinitely slow evolution, namely, adiabatically. Equation (3) includes also the limit 
t — > — oo. The order of the two limits is of course essential and cannot be interchanged. 

Intuitively the content of the Gell-Mann and Low theorem is simple: if the 
Hamiltonian evolves adiabatically and if we start from the ground state of the 
Hamiltonian H(to) at a given initial time to, the system will remain in the ground 
state of the local Hamiltonian H(t) at any subsequent time t, since an infinitely slow 
evolution cannot excite any system by a finite amount of energy. It is therefore essential 
for the validity of the theorem that, during the evolution, the local ground state never 
becomes degenerate, e.g. no phase transition occurs. In the latter case, Eq. (3) will 
provide a state ipo which is not the ground state of H but the state which can be obtained 
smoothly from the unperturbed ground state O through the adiabatic switching on of 
the interaction. 

The operator U^(t) can be obtained by a perturbation expansion from the free 
evolution operator Uo(t) = exp(—iH t/h), and for the present purpose one can write 

U(e) = 1 - i J ^ Hj^dh + {-If J ^ Hj(t 2 )dt 2 Hjitjdh ■ ■ ■ 

/-oo dhT [if 7 (t n )#j(*„-i) • • • #/(ti)] (5) 

where T is the time ordered operator and 

Hj(t) = e iHot/h H e^ e -iH Q t/h _ (q) 

In Eq. (6) the indication of the dependence of Hj on e was omitted for simplicity. The 
limit e — > has to be taken after all the necessary manipulations have been performed. 
The demonstration of the Gell-Mann and Low theorem, based on the expansion of Eq. 
(5), can be found in the original paper or in textbooks on general many-body theory 
(Fetter and Walecka 1971). From Eq. (3), it follows that the energy shift AE due to 
the nucleon-nucleon interaction is given by 

AE = u m <t'z ( ;' ( ~°^ o> . p> 

6-o (0 o |^)(-oo)|0o) 
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where the expansion of Eq. (5) has to be used both in the numerator and in the 
denominator. The procedure is ill-defined in the limit e — > 0, as one can see by 
considering the first non-trivial term (n = 1) of the expansion of Eq. (5) and taking 
the matrix elements appearing in Eq. (7). They blow up in that limit. Fortunately, 
here we can get help from the so called "linked cluster" theorem. The formulation of 
the theorem is better stated in the language of the diagrammatic method, as explained 
below. The theorem shows that the numerator and the denominator possess a common 
factor, which includes all the diverging terms, and therefore they cancel out exactly, 
leaving a well defined finite result. 

Finally, each term of the perturbation expansion can be explicitly worked out by 
means of Wick' s theorem, which allows one to evaluate the mean value of an arbitrary 
product of annihilation and creation operators in the unperturbed ground state. Then 
the perturbative expansion of the interaction energy part AE of the ground state energy 
can be expressed in terms of "Goldstone diagrams", as devised by Goldstone (1957). 
Each diagram represents, in a convenient graphical form, a term of the expansion, in 
order to avoid lengthy analytical expressions and to make their structure immediately 
apparent. The general rules (from i to vi below) for associating the analytical expression 
to a given diagram are described in the following. The expression is constructed by the 
following factors. 

(i) Each drawing of the form of Fig. 1, which can be called conventionally a "vertex", 
as usual represents a matrix element of the two-body interaction, according to the rules 
discussed previously. 

(ii) A line with an upward (downward) arrow indicates a particle (hole) state, and it 
will be labeled by a momentum k (including spin-isospin), a different one for each line. 

(in) Between two successive vertices a certain number of lines (holes or particles) will 
be present in the diagrams. Then, the energy denominator 



e Efc, E k . - J2 k [ E k[ +*V 

is introduced, where now the summation runs only on the particle and hole energies 
which are present in the diagram between the two vertices. 

(iv) Each diagram is given an overall sign (-i)'H-H-n-i^ where n is the order of the 
diagram in the expansion, h is the total number of hole lines in the diagram, and I the 
number of closed loops. A "loop" is a fermion line (hole or particle) which closes on 
itself when followed along the diagrams, as indicated by the directions of the arrows, 
passing eventually through the dots of vertices. 

(v) Finally a "symmetry factor" of the form (|) s , s = 0, 1, 2 • • •, has to be put in front 
of the whole expression. In general, the factor is connected with the symmetry of the 
diagram, and to find its correct value it is necessary to analyze the formalism in more 
detail. Let us consider the case where two lines, both particles or both holes, connect the 
same two interaction vertices, without being involved in any other part of the diagram. 
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They can be called "equivalent lines" . In this case, the only one that will be considered, 
the symmetry factor is \. 

(vi) Of course, one must finally sum over all the momenta labeling the lines of the 
diagram. 

Since we are considering the ground state energy, only "closed" diagrams must be 
included, i.e. no external line must be present. Furthermore, according to the linked- 
cluster theorem, only connected diagrams must be considered, i.e. diagrams which 
cannot be separated into two or more pieces with no line joining them. 

In conclusion, the ground state energy shift is obtained by summing up all possible 
closed and linked diagrams 

AE = lim(0 o |/7 1 ^(-oo)|0o)cL , (9) 

where the subscript CL means connected diagrams only, linked with the first interaction 
Hi. The latter specification means simply, in this case, that the diagram must be 
complete, namely it must involve all the interactions. 

Another fundamental consequence of the restriction to connected diagrams is that 
the energy shift AE is proportional to the volume of the system, as it must be for 
extended systems with short range interactions only. Disconnected diagrams have the 
unphysical property to be proportional to higher powers of the volume. 

Let us consider the nuclear matter case with a typical NN interaction. The NN 
interaction is characterized by a strong repulsion at short distance. The simplest 
assumption would be to consider an infinite hard core below a certain core radius. Such 
a NN potential has obviously infinite matrix elements in momentum representation, 
and a perturbation expansion has no meaning. All modern realistic NN interactions 
introduce a finite repulsive core, which however is quite large, and therefore in any 
case a straightforward perturbative expansion cannot be applied. The repulsive core is 
expected to modify strongly the ground state wave function whenever the coordinates of 
two particles approach each other at a separation distance smaller than the core radius 
c. In such a situation the wave function should be sharply decreasing with the two 
particle distance. The "wave function" of two particles in the unperturbed ground state 
0o can be defined as (ki, k 2 < hp) 

0(n,r 2 ) = (001^1)^2)^^00) = e^+^M^W 2 , (10) 

where (i are spin-isospin variables, and R = (ri + r 2 )/2, r = (ri — r 2 ) are the center 
of mass and relative coordinate of the two particles respectively. Therefore the wave 
function of the relative motion in the s-wave is proportional to the spherical Bessel 
function of order zero jo(kr), with k the modulus of the relative momentum vector 
k = (ki — k 2 )/2. The core repulsion is expected to act mainly in the s-wave, since it 
is short range, and therefore this behaviour must be strongly modified. In the simple 
case of k = the free wave function jo(kr) — > 1, and schematically one can expect a 
modification, due to the core, as depicted in Fig. 2. The main effect of the core is to 
"deplete" the wave function close to r = 0, in a region of the order of the core radius 
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c. Of course, the attractive part of the interaction will modify this simple picture at 
r > c. If the core interaction is the strongest one, then the average probability p for 
two particles to be at distance r < c would be a measure of the overall strength of the 
interaction. If p is small, then one can try to expand the total energy shift AE in power 
of p. The power p n has, in fact, the meaning of probability for n particles to be all at 
a relative distance less than c. In a very rough estimate p is given by the ratio between 
the volume occupied by the core and the average available volume per particle 

P « (s) 3 in) 

with ^-d 3 = p^ 1 . From Eq. (11) one gets p « ^(/cpc) 3 , which is small at saturation, 
= 1.36 /to -1 , and the commonly adopted value for the core is c = 0.4 /to -1 . The 
parameter remains small up to few times the saturation density. 

The graphs of the expansion can now be ordered according to the order of the 
correlations they describe, i.e. the power in p they are associated with. It is easy 
to recognize that this is physically equivalent to grouping the diagrams according 
to the number of hole lines they contain, where n hole lines correspond to n-body 
correlations. In fact, an irreducible diagram with n hole lines describes a process in 
which n particles are excited from the Fermi sea and scatter in some way above the 
Fermi sea. Equivalently, all the diagrams with n hole lines describe the effect of clusters 
of n particles, and therefore the arrangement of the expansion for increasing number of 
hole lines is called alternatively "hole expansion" or "cluster expansion". 

The series of two hole-line diagrams starts with the diagrams depicted in Fig. 3 
(first order) and in Fig. 4 (second order) and continues with the ones shown in Fig. 
5. The infinite set of diagrams depicted in Fig. 5 can be summed up formally by 
introducing the two-body scattering matrix G, as schematically indicated in Fig. 6. In 

A 



0(r) 



c 

r 

Figure 2. Schematic representation of the expected effect of the core repulsion on the 
two-body wave function in nuclear matter. 




Figure 3. Direct and exchange first order diagrams. 
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Figure 6. The geometric series for the G-matrix. 



Equation of State of Nuclear Matter at high baryon density 



9 




k 2 

Figure 7. Two hole- line diagrams for the ground state energy in terms of the G-matrix 
(wiggly lines). 

the second line of Fig. 6 the geometric series has been re-introduced, once the initial 
interaction has been isolated. This corresponds to the following integral equation 

(kik 2 \G(u)\k s k4) = (hk 2 \v\k 3 k 4 ) + 

+ E KK (kik2\v\k' 3 k' 4 ) (l " e ^ e y^ FW)) (m\G(u)\k 3 k 4 ) . (12) 

In the diagrams, the intermediate states are particle states, and this is indicated in 
Eq. (12) by the two factors 1 — Qp(k). One can consider the diagrams of Fig. 6 
as part of a given complete diagram of the total energy expansion. Therefore all the 
energy denominators contain the otherwise undefined quantity iv, usually indicated as 
the "entry energy" of the G matrix. The precise value of u> will depend on the rest 
of the diagram where the G matrix appears, as we will see soon. Equation (12) is 
anyhow well defined for any given value of u>. It has to be noticed that Eq. (12) is very 
similar to the equation which defines the usual off-shell scattering T matrix between two 
particles in free space. The G matrix of Eq. (12) can be considered the generalization 
of the T matrix to the case of two particles in a medium (nuclear matter in our case). 
Actually in the zero density limit, Qp(k) — > and the G matrix indeed coincides with 
the scattering T matrix. Once the G matrix has been introduced, the full set of two 
hole-line diagrams can be expressed as in Fig. 7, where the G matrix is now indicated 
by a wiggly line. This notation stresses the similarity between the G matrix and the 
bare nucleon-nucleon interaction v. This result, as depicted in Fig. 7, can be checked 
by expanding Eq. (12) (by iteration). The entry energy in this case is u> = + ek 2 , 
which means that the G matrix is "on the energy shell" , i.e. the G matrix is calculated 
at the energy of the initial state. The diagrams need a factor |, since the two hole-lines 
are equivalent, according to rule (v). Therefore the correction AE 2 to the unperturbed 
total energy (just the kinetic energy), at the two hole-line level of approximation, is 
given by 

A£ 2 = l - Y, (hk 2 \G(e kl + e k2 )\hk 2 } A , (13) 

ki,k'2<kp 

where the label A indicates that both direct and "exchange" matrix elements have 
to be considered, i.e. \kik 2 )A = \kik 2 ) — \k 2 ki). One of the major virtues of the G matrix 
is to be defined even when the interaction v is singular (e.g. it presents an infinite hard 
core). This shows that the G matrix is in some sense "smaller" than the NN interaction 
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v, and an expansion of the total energy shift AE in G, instead of v, should have a better 
degree of convergence. To substitute v with the matrix G in the original expansion is 
always possible, since a "ladder sum" (a set of diagrams of the type in Fig. 6) can 
always be inserted at a given vertex and the corresponding series of diagrams summed 
up (with the proviso of avoiding double counting). In general, however, the resulting 
G-matrix will be "off the energy shell" , which complicates the calculations considerably. 
It turns out, anyhow, that also the bare expansion of AE in terms of the G-matrix, in 
place of the NN interaction v, is still badly divergent. 

The solution of this problem is provided by the introduction of an "auxiliary" single 
particle potential U (k). The physical reason of such a procedure becomes apparent if one 
notices that the energies of the hole or particle states are surely modified by the presence 
of the interaction Hi, and intuitively they should have some relevant effects on the total 
energy of the system. However, in the Goldstone expansion of Eq. (9), or similar, such 
an effect does not appear explicitly, and therefore it should be somehow introduced into 
(or extracted from) the expansion, since, physically speaking, any two-body or higher 
correlations should be evaluated as corrections to some mean field contribution. The 
genuine strength of the correlations has to be estimated in comparison with a reference 
mean field energy, rather than to the free particle energy. The explicit form of the 
auxiliary single particle potential has to be chosen in such a way to minimize the effect 
of correlations, which is equivalent to speed up the rate of convergence of the expansion. 
Formally, one can re-write the original Hamiltonian by adding and subtracting the 
auxiliary single particle potential U 

H = (H + U) + (Hi -U)=H' Q + H[ 

= Hk^ka\a k , (14) 

and consider e k as the new single particle spectrum. The expansion is now in the 
new perturbation interaction H[. The final result should be, of course, not dependent 
on U, at least in principle. A "good" choice of the auxiliary potential U is surely 
one which is able to strongly reduce the contribution of H[ to the total energy of the 
system. The perturbation expansion in H[ can be formulated in terms of the same 
Goldstone diagrams discussed previously, where the single particle kinetic energies t k 
are substituted by the energies e k = t k + U (k) in all energy denominators. Furthermore, 
new terms must be introduced, which correspond to the so called U U insertions". More 
precisely, the rules i-vi above must be supplemented by the following two other additional 
rules. 

(i - bis) A symbol of the form reported in Fig. 8 indicates a U insertion, which 
corresponds to a factor [7(/ci)5fc(ki — k 2 )%£ 2 in the diagram. 

(iv - bis) A diagram with a number u of U insertions contains the additional phase 
(— l) u . This is a trivial consequence of the minus sign with which U appears in H[. 

The first U insertion of Fig. 9 cancels out exactly the potential energy part of 
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Figure 8. Representation of a potential insertion factor. 




Figure 9. The first potential insertion diagram. 

the single particle energies as contained in H' , see Eq. (14), and therefore the total 
energy at the two hole-line level is given by 

E ^=T,^ + l E (kMG^+e^lhk^A . (15) 

k<kp ki,k2<kp 

One has to keep in mind that the G matrix depends now on U, since the auxiliary 
potential appears in the definition of the single particle energies e^. The appearance 
of the unperturbed kinetic energy is valid for any choice of the auxiliary potential and 
it is not modified by the addition of the higher order terms in the expansion. It is 
a distinctive feature of the Goldstone expansion that all correlations modify only the 
interaction part and leave the kinetic energy unchanged. Of course this property is 
pertinent only to the expression of the ground state energy. 

It is time now to discuss the choice of the auxiliary potential. A good choice of U 
should minimize the contributions from higher order correlations, i.e. the contributions 
of the diagrams with three or more hole-lines. In other words, the U insertion diagrams 
must counterbalance the diagrams with no U insertion. An exact cancellation is of course 
not possible, however one can select some graphs which are expected to be large and try 
to cancel them out exactly. At the three hole-line level, one of the largest contributions 
is expected to be given by the graph of Fig. 10a. In this diagram the symbol, already 
introduced, for indicating a G matrix stands for the corresponding ladder summation 
inside the diagram. This can be done systematically along the expansion, but one has 
to be careful in checking the energy at which the G matrix has to be calculated, i.e. if 
it is on shell or off shell. It has been shown by Bethe, Brandow and Petschek (1962) 
that it is possible to choose U in such a way that the corresponding potential insertion 
diagram, shown in Fig. 10b, cancels out the (hole) "bubble diagram" of Fig. 10a. This is 
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Figure 10. One of the lowest order (in the G-matrix) three hole-line diagrams (a) 
and the corresponding potential insertion diagram (b). 

indeed possible by virtue of the so called BBP theorem established by the authors, which 
states that the G matrix connected with the bubble in the diagram of Fig. (10a) must 
be calculated on the energy shell, namely u = e kl + e k2 . For the other two G matrices 
appearing in the diagram this property is also valid, but this is a trivial consequence of 
the theorem. Therefore, if one adopts for the auxiliary potential the choice 

U(k) = J2 (kk'\G(e kl +e k2 )\kk') , (16) 

k'<k F 

it is straightforward to see that the diagram of Fig. 10b is equal to minus the diagram 
of Fig. 10a (remind the rule iv - bis). 

The choice of Eq. (16) for U was originally devised by Brueckner, on the basis 
of physical considerations. The choice of Eq. (16) is therefore called the Brueckner 
potential; it implies a self-consistent determination of U, since, as already mentioned, 
the G matrix itself depends on U. The hole expansion with the Brueckner choice for U 
is called the Bethe-Brueckner-Goldstone (BBG) expansion. 

In the original Brueckner theory the potential U was assumed to be zero above hp. 
This is called the "standard choice", or "gap choice", since it necessarily implies that 
the single particle energy e k is discontinuous at k — kp. This choice also implies that 
the potential insertion diagram of Fig. lib is automatically zero. The corresponding 
diagram, with the G matrix replacing the auxiliary potential, depicted in Fig. 11a, is 
therefore in no way counterbalanced. The G matrix in this diagram is off shell. In fact, 
the BBP theorem does not hold for it. The graph of Fig. 11a is usually referred to also 
as the particle bubble diagram, or simply "bubble diagram", and in the following this 
terminology is adopted. 

Another possible choice for the auxiliary potential U(k) is the so called "continuous 
choice", where U(k) is defined by Eq. (16) for all values of |k|. In this case the 
potential is continuous through the Fermi surface and e(k) can be interpreted as a 
single particle spectrum. Furthermore the two diagrams of Fig. 11 can have some 
degree of compensation, as we will see in the applications. Since the final results must 
be independent of the choice of the auxiliary potential, the sensitivity of the results 
to U{k) at a given order of the expansion can be used as a criterion for the degree of 
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Figure 11. Particle bubble diagram (a) and the corresponding potential insertion 
diagram (b). 



convergence reached at that level of approximation. No sensitivity would correspond to 
a complete convergence. 

The bubble diagram of Fig. (11a) can be considered the first term of the full 
set of three hole-line diagrams. The two hole-line diagrams have been summed up by 
introducing the two-body G matrix, which is the generalization to the nuclear medium 
of the two-body scattering matrix in free space. From Eq. (12) it is apparent that 
the only difference between the G matrix and the free space scattering matrix is the 
presence of the "Pauli operator" Q(k 1: k 2 ) = (1 — 0f(/ci))(1 — Q F (k 2 )), with 6^(/c) the 
(zero temperature) Fermi distribution, and the presence of the energies e& in place of 
the kinetic energies. This has far-reaching consequences. 

It is therefore conceivable that the three hole-line diagrams could be summed up 
by introducing some similar generalization of the scattering matrix for three particles 
in free space, which would correspond physically to consider the contribution of the 
three-body clusters. The three-body scattering problem has a long history by itself, and 
has been given a formal solution by Fadeev (1965). For three distinguishable particles 
the three-body scattering matrix is expressed as the sum of three other scattering 
matrices, = T 1 +T 2 +T 3 . The scattering matrices T; satisfy a system of three coupled 
integral equations. The kernel of this set of integral equations contains explicitly the 
two-body scattering matrices pertaining to each possible pair of particles. Also in this 
case, therefore, the original two-particle interaction disappears from the equations in 
favor of the two-body scattering matrix. The formal reason for this substitution is the 
need of avoiding "disconnected processes" , which introduce spurious singularities in the 
equations (Fadeev 1965). For identical particles the three integral equations reduce to 
one, because of symmetry. In fact, the three functions Tj must coincide within a change 
of variable with a unique function, which we can still call T^ 3 \ The analogous equation 
and scattering matrix in the case of nuclear matter (or other many-body systems in 
general) has been introduced by Rajaraman and Bethe (1967). The integral equation, 
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Figure 12. Expansion of the Bethe-Fadeev integral equation. 



the Bethe-Fadeev equation, reads schematically 

T^ 3 ) = G + G X Ql y( 3 ) 

e 

(k 1 k 2 k 3 \T^\k' 1 k' 2 k' 3 ) = (k 1 k 2 \G\k' 1 k' 2 )5 K (k 3 -k' 3 ) + 

+ (k 1 k 2 h\G 12 X^T^\kikf 2 k' 3 ) . (17) 

The factor Q 3 /e is the analogous of the similar factor appearing in the integral equation 
for the two-body scattering matrix G, see Eq. (12). Therefore, the projection operator 
Q 3 imposes that all the three particle states lie above the Fermi energy, and the 
denominator e is the appropriate energy denominator, namely the energy of the three- 
particle intermediate state minus the entry energy u, in close analogy with the equation 
for the two-body scattering matrix G, Eq. (12). The real novelty with respect to the 
two-body case is the operator X . This operator interchanges particle 3 with particle 
1 and with particle 2, X = Pi 23 + P± 32 , where P indicates the operation of cyclic 
permutation of its indices. It gives rise to the so-called "endemic factor" in the Fadeev 
equations, since it is an unavoidable complication intrinsic to the three-body problem 
in general. The reason for the appearance of the operator X in this context is that no 
two successive G matrices can be present in the same pair of particle lines, since the 
G matrix already sums up all the two-body ladder processes. In other words, the G 
matrices must alternate from one pair of particle lines to another, in all possible ways, as 
it is indeed apparent from the expansion by iteration of Eq. (17), which is represented 
in Fig. 12. 

Therefore, both cyclic operations are necessary in order to include all possible processes. 
In the structure of Eq. (17) the third particle, with initial momentum k 3 , is somehow 
singled out from the other two. This choice is arbitrary, but it is done in view of the 
use of the Bethe-Fadeev equation within the BBG expansion. 

In order to see how the introduction of the three-body scattering matrix allows 
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Figure 13. Schematic representation of the direct (a) and exchange (b) three hole-line 
diagrams. 

one to sum up the three hole line diagrams, we first notice, following Day (1981), that 
this set of diagrams can be divided into two distinct groups. The first one includes the 
graphs where two hole lines, out of three, originate at the first interaction of the graph 
and terminate at the last one without any further interaction in between. Schematically 
the sum of this group of diagrams can be represented as in Fig. 13a. The third hole line 
has been explicitly indicated, out from the rest of the diagram. The remaining part of 
the diagram describes the scattering, in all possible ways, of three particle lines, since no 
further hole line must be present in the diagram. This part of the diagram is indeed the 
three-body scattering matrix T^ 3 \ and the operator Q 3 in Eq. (17) ensures, as already 
mentioned, that only particle lines are included. 

The second group includes the diagrams where two of the hole lines enter their 
second interaction at two different vertices in the diagram, as represented in Fig. 13b. 
Again the remaining part of the diagram is T&\ i.e. the sum of the amplitudes for all 
possible scattering process of three particles. It is easily seen that no other structure 
is possible. The set of diagrams indicated in Fig. 13b can be obtained by the ones of 
Fig. 13a by simply interchanging the final (or initial) point of one of the "undisturbed" 
hole lines with the final (or initial) point of the third hole line. This means that one 
can obtain each graph of the group depicted in Fig. 13b by acting with the operator 
X on the bottom of the corresponding graph of Fig. 13a. In this sense the diagrams 
of Fig. 13b can be considered the "exchange" diagrams of the ones in Fig. 13a (not to 
be confused with the term "exchange" previously introduced for the matrix elements 
of G). If one inserts the terms obtained by iterating Eq. (17) inside these diagrams in 
substitution of the scattering matrix (the box in Fig. 13), the first diagram, coming 
from the inhomogeneous term in Eq. (17) is just the bubble diagram of Fig. 11a. The 
corresponding exchange diagram is the so called "ring diagram" , reported in Fig. 14. It 
turns out that for numerical reasons it is convenient to separate both bubble and ring 
diagrams from the rest of the three hole-line diagrams, which will be conventionally 
indicated as "higher" diagrams. 

Indeed, going on with the iterations, one gets sets of diagrams as the ones depicted 
in Figs. 15, and so on. 



Equation of State of Nuclear Matter at high baryon density 16 




Figure 14. The ring diagram, belonging the set of three hole-line diagrams. 




Figure 15. The series of three hole-line diagrams, up to fifth order in the G-matrix. 

To these series of diagrams one has, of course, to add the diagrams obtained by 
introducing the exchange matrix elements of G in place of the direct ones (if they really 
introduce a new diagram). The structure of the diagrams of Figs. (15) displays indeed 
the successive three-particle scattering processes. 

Let us notice that the graph of Fig. 10a, where the bubble is attached to the 
hole line, is not included, and it has to be added separately, as previously discussed in 
connection with the U insertion diagrams. 

Some ambiguity arise if the diagram of Fig. 16 should be included at the three hole- 
line level or not. The diagram is usually referred to as the "hole-hole" diagram, for 
obvious reasons. Although, due to momentum conservation, only three hole lines are 
independent, we will consider this particular diagram as belonging to the four hole-line 
class. 

For writing down explicitly the three hole-line contribution to the total energy we 
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Figure 16. The hole-hole diagram. It is not included in the Bethe-Fadeev equation. 



+ 




Figure 17. Direct and exchange contribution within the three hole-line diagrams. 



still need to find out the correct symmetry factors and signs. Let us first consider the 
part of the diagram which describes the interaction among the three particle lines. In 
the scattering processes each two-body G matrix can involve both the direct and the 
exchange term, as illustrated in Fig. 17. Hence, there is no additional symmetry factor 
involved. The three lines, in fact, are never equivalent, since the various G matrices are 
alternating among the different possible pairs of particles along the diagram. Therefore, 
the direct and exchange matrix elements of each G matrix have to be considered and 
no symmetry factor for this part of the diagram has to be introduced. Let us consider 
now the hole lines which close the diagram. For the diagrams of the type of Fig. 13a, 
two equivalent hole lines appear, joining the first and the last interaction. As discussed 
previously, this implies the introduction of a symmetry factor equal to \ in front of 
each diagram belonging to this group. In conclusion, the explicit expression for the 
contribution of the whole set of diagrams of Fig. (13a) (the "direct" diagrams) can be 
written 

^th = lJ2k 1 toM<kF^{k'},{k"}>k F {klk2\G\k' 1 k' 2 )A 

■\ {k' x k' 2 k^\XT^X\k'[k'{k'i) ± {k'[k'{\G\k x k 2 ) A , (18) 

where again the operators X are introduced in order to generate all possible scattering 
processes, with the condition that the G matrices alternate, from one interaction to 
the next one, among the possible pairs out of the three particle lines. In Eq. (18) the 
denominator e = E k i + E k i — E kl — Ek 2 , and analogously e' = E k rr + E k » — E kl — E k , 2 . 
Let us now consider the exchange diagrams of Fig. 13b. They can be obtained by 
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interchanging the initial point (or end point) of the hole line labeled k% with the one 
of the hole line k 2 , i.e. by multiplying by P 12 3 the expression of Eq. (18) (at the right 
or left side). In this case, however, we have to omit the symmetry factor |, since no 
pair of equivalent lines appears any more. We could equally well interchange k 3 with 
ki, since in this way we actually consider the same set of diagrams. This can be readily 
checked by displaying explicitly the sets of associated diagrams in the two cases. It 
is then convenient to take the average of the two possibilities, which is equivalent to 
multiply by P 12 3 + P132 = X the expression of Eq. (18) and to reintroduce the factor \. 
In summary, the entire set of three hole-line diagrams can be obtained by multiplying 
the expression of Eq. (18) by 1 + X. 

It is convenient in Eqs. (17) and (18) to single out the first interaction which occurs 
in T^ 3 \ where the third hole line must originate. Posing = GD, or, explicitly 

(hhhlT^lk'^k's) = ^( k MG\k'(k^(k';k 2 %\D\k[k 2 k' 3 ) , (19) 

h" h" 

then the matrix D satisfies the formal equation 

D — 1 — X—GD . (20) 

e 3 

Notice that, contrary to the G-matrices appearing in Eq. (13), the G matrix appearing 
in Eq. (20) is off-the energy shell, since the denominators which enter in its definition 
contain the energy of the hole lines k 1: k 2) k 3 , as well as of the third particle line. The 
denominator e 3 is the energy of the appropriate three particles-three holes intermediate 
state. 

Summarizing, the three-hole line contribution can be obtained by solving the 
integral equation (20) and inserting the solution in Eq. (18). Notice that the solution 
D depends parametrically on the external three hole lines momenta. 

The scattering matrix (or equivalently D) can be used as the building block 
for the construction of the irreducible four-body scattering matrix in an analogous 
way as the two-body scattering matrix G has been used to construct T^ 3 \ The resulting 
equations for in the case of four particles in free space, are called Yakubovsky 
equations. It is not difficult to imagine that additional complexities are involved in 
these equations. Since nobody till now has dared to write down these equations for 
nuclear matter, not to say to solve them, we will not discuss their structure. However, 
estimates of the four-hole lines contribution have been considered (Day 1981). 

3. Nuclear matter within the BBG expansion. 

Before summarizing the theoretical results for the EoS at zero temperature on the 
basis of the BBG expansion, let us briefly analyze in more detail the properties of the 
scattering matrix G. As already mentioned, the G matrix can be considered as the in 
medium two-body scattering matrix. This can be more clearly seen by introducing the- 
two body scattering wave function ^k lt k 2 in analogy to the case of free space scattering 
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(Newton 1966) 

|*fc 1)fc2 ) = \kifa) + -G\hk 2 ) = \hk 2 ) + -v\* klM ) , (21) 
" e e 

where we have used the relationship G\k\k 2 ) = v^^m)- The latter is obtained by 

multiplying by v the first of Eqs. (21), which defines the scattering wave function |\&) , 

and making use of the integral equation (12) for the G matrix. It is instructive to look 

more closely to the scattering wave function in coordinate representation. The centre 

of mass motion separates, since the total momentum P is a constant of the motion. In 

the notation of Eq. (10), the integral equation for the scattering wave function, in the 

relative coordinate, reads 

^(r) = e ikr + / dV(r|^|r>(r)V>(r') , (22) 
J e 

where, for simplicity, the spin-isospin indices have been suppressed and the NN 
interaction has been assumed to be local. Still the wave function ip depends on both the 
total momentum P and on the entry energy u. The latter appears in the denominator 
e, see Eq. (12), while the total momentum P appears also in the Pauli operator Q. The 
kernel ^ in Eq. (22) is the same as in the usual theory of two-body scattering (Newton 
1966), except for the Pauli operator Q, which has a deep consequence on the properties 
of ip. This can be most easily seen if one considers the case P = and an entry energy 
corresponding to two particles inside the Fermi sphere, ui < 2Ep. In this case the Pauli 
operator simply implies that the relative momentum |k|' > hp, and the kernel reads, 
after a little algebra 

^i^i^- 1 r°° k ' dk ' smfe/ i r - r 1 

(T| — \T ) - — / — - - — , [Z6) 



e 2ir 2 Jk F ley — u 

where the energy denominator never vanishes (provided the energy ey is an increasing 
function of k', as it always happens in practice). In the usual scattering theory (Newton 
1966), on the contrary, the denominator can vanish and the integral on k' provides the 
free one particle Green' s function (according to the chosen boundary conditions). Then, 
for large values of r, one gets the usual asymptotic behaviour (for outgoing boundary 
condition) 

V(r) - e^-W)— , (24) 
r 

which describes an outgoing spherical wave and therefore a flux of scattered particles. 
Here 9 is the angle between r and the initial momentum k. The asymptotic behaviour 
of the kernel of Eq. (22) can be obtained by a first partial integration with respect to 
the sine function 

v 1 e 1 ; 2n 2 2e kF - u |r-r'| 2 Vl y ' 

since further partial integrations give higher inverse power of |r — r'| (here the non- 
vanishing of the energy denominator is essential). The asymptotic behaviour of ip(r) 
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follows easily from Eq. (25), since the NN interaction is of short range. Inserting Eq. 
(25) in Eq. (22), one gets 

cos /c^r 



ikr 



V>(r) - e lk - 



(26) 



In this case the scattered flux vanishes at large distance, since the scattered wave 
vanishes faster than 1/r, and no real scattering actually occurs. The scattering wave 
function ^(r) indeed merges, at large distance r, into the two-body relative wave function 
of Eq. (10) for a gas of free particles. In the language of scattering theory this means 
that all the phase shifts are zero. This property is usually called the "re-phasing" of 
the function ip. The two-body wave function does not describe a scattering process 
but rather the distortion of the two-body relative motion due to the interaction with 
respect to free gas case. Since the interaction is assumed to be of short range, such 
a distortion is concentrated at short distance, mainly inside the repulsive core region, 
but also slightly outside it (due to the attractive part of the interaction and to quantal 
effects). 

It has to be stressed that for entry energy corresponding to two particles above 
the Fermi sphere, the two-body wave function ip(r) can be still defined, as well as 
the scattering G matrix, and this is indeed necessary for the continuous choice of the 
auxiliary potential U. In this case ip(r) can describe a real scattering (the energy 
denominator can vanish), namely a collision process of two particles inside nuclear 
matter, provided the correct boundary condition is imposed. This is always the case 
when the two particles initial momenta lie above the Fermi sphere. 
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Figure 18. Two-body relative wave-function for the free nucleon gas (full line) and 
for correlated nuclear matter (dots). 



Let us go back to the case of two particles inside the Fermi sphere. The difference 
between the wave function ip(r) and the corresponding free wave function exp(i(k • r)), 
already introduced, Eq. (10), is called the "defect function" Cfei,fc 2 - The size of the 
defect function is a measure of the two-body correlations present in the system. As a 
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more quantitative parameter one can take the norm of the defect function, averaged 
over the Fermi sphere and calculated inside the available volume per particle. The 
parameter is usually called "wound parameter", since it describes the "wound" in the 
wave function produced by the NN correlations, and it is a more refined version of 
the parameter p previously introduced in discussing the hole expansion. Since the 
repulsive core is expected to have the dominant effect in nuclear matter, and it is of 
short range, the s-wave component of ip(r) should be the most affected one, and therefore 
the corresponding defect function should be the largest one. This is indeed the 
shown in Fig. 18, where the function ip(r) in the channel 1 S , (Eqs. (21) and (22) can be 
easily written in spin-isospin coupled representation), is reported (dots) in comparison 
with the corresponding free relative wave function (full line). The calculations have 
been done in the continuous choice and at saturation density &f = 1.36 /to -1 . The 
initial relative momentum was chosen at q — 0.1 /to -1 (the value at q = exactly can 
create numerical problems). One can see that the distortion of the free wave function 
is concentrated at small r values (the square of the wave function has to be taken). 
One can notice the striking similarity with the naive guess of Fig. 2. At distance 
larger than the core radius the correlated wave function oscillates slightly around the 
uncorrelated one, an effect mostly due to the large distance attractive component of the 
NN interaction. 

The short distance correlation is expected to decrease for higher partial waves. It 
should also be affected by the initial relative momentum. Both effects are shown in Fig. 
19, where the two-body wave function at relative momentum q = &p is reported for the 
1 S , 1 P\ and 1 D 2 channels. The "healing effect" is apparent in all these cases, the two- 
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Figure 19. The same as in Fig. 18, but at relative momentum q — kp and for the 
three channels x 5o, 1 P\ and 1 D 2 - The first peak of the wave function is decreasing at 
increasing values of the partial wave I = 0,1,2. 



body relative wave function is strongly suppressed at short distance. The corresponding 
"healing distance" , the size of the interval where the suppression occurs, can be used 
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to estimate the wound parameter, as discussed above. Values of this parameter are 
about 0.2-0.25 around saturation density for symmetric nuclear matter, which indicate 
a moderate rate of convergence. Even at densities of few times the saturation value 
the wound parameter does not exceed 0.3 - 0.35. In pure neutron matter it turns out 
that the wound parameter is smaller by about a factor 2 in the same density range, and 
convergence should be much better. 

It has to be stressed, anyhow, that the reduction of weight that should be obtained 
by an additional hole line in the set of diagrams along the BBG expansion does not 
depend only on the probability to find two particles at short distance, but also on the 
action of the NN potential on the defect function (. Since the introduction of the 
scattering G-matrix should take care, to a large extent, of the short range correlations 
due to the repulsive core, the higher order correlations should contain a more balanced 
contribution from the repulsive and attractive parts of the interaction, and therefore a 
strong compensation between attractive and repulsive contributions to the expansion. 
With some degree of optimism, one can hope that the expansion rate could be even 
better than the one guessed from the value of the wound parameter. 

This expectation is indeed confirmed by actual calculations of the three hole-line 
contributions. The results of Baldo et al. (2001) are reported in Fig. 20 for the Argonne 
Vis NN potential (Wiringa et al. 1995), and symmetric nuclear matter, both for the 
gap and for the continuous choice. The full lines correspond to the Brueckner two 
hole-line level of approximation (Brueckner-Hartree-Fock or BHF), while the symbols 
indicate results obtained adding the three hole-line contributions. Two conclusions can 
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Figure 20. Equation of state of symmetric nuclear matter at the two hole-line level 
(full lines) in the gap (BHF-G) and in the continuous choice (BHF-C) of the single 
particle potential. The symbols label the corresponding EoS when the three hole-line 
contributions are added. 



be drawn from these results. 
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i) At the Brueckner level the gap and continuous choice still differ by few MeV, which 
indicates that the results depend to a certain extent on the choice of the auxiliary 
potential. According to the discussion above this implies that the expansion has not 
yet reached full convergence. On the contrary when the three hole-line diagrams are 
added the results with the different choices for the single particle potential U are quite 
close, which is surely an indication that the expansion has reached a good degree of 
convergence. Notice that the insensitivity to the choice of U is valid in a wide range of 
density, only at the highest density some discrepancy starts to appear. One can see that 
even at 4-5 times saturation density the BBG expansion can be considered reliable. 

ii) As already discussed, the auxiliary potential is crucial for the convergence of the 
BBG expansion. It is not surprising, therefore, that the rate of convergence is dependent 
on the particular choice of U . From the results it appears that the continuous choice is 
an optimal one, since the three hole-line corrections are much smaller and negligible in 
first approximation. 

It is important to stress that the smallness of the three hole-line corrections is the 
result of a strong cancellation of the contributions of the different diagrams discussed 
above. This is illustrated in Fig. 21, where the values of the bubble (figure 11a), ring 
(figure 14), U-potential insertion (lib) and higher order diagrams are reported. This 
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Figure 21. The contributions of different three hole-line diagrams to the interaction 
energy in symmetric nuclear matter. The curve "bub-U" gives the contribution of the 
potential insertion diagram of Fig. lib. The line labeled "total" is the sum of all 
contributions. 

shows clearly the relevance of grouping the diagrams according to the number of hole 
lines, in agreement with the BBG expansion. An ordering of the diagrams according to 
e.g. the number of G-matrices involved would be badly divergent. 

Similar results are obtained for pure neutron matter, as illustrated in Fig. 22, taken 
from Baldo et al. (2000). Notice that, in agreement with the previous discussion on the 
wound parameter, the rate of convergence looks faster in this case. 
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Figure 22. Equation of state of pure neutron matter at the two hole-line level 
(full lines) in the gap (BHFG) and in the continuous choice (BHFC) of the single 
particle potential. The symbols label the corresponding EoS when the three hole-line 
contributions are added. 



4. The Coupled Cluster Method 

The BBG expansion can be obtained also within the Coupled-Cluster Method (CCM) 
(Kummel et al. 1978), a general many-body theory which has been extensively applied 
in a wide variety of different physical systems, both bosonic and fermionic ones. The 
connection between the CCM and the BBG expansions has been clarified by Day (1983). 
In the CCM method one starts from a particular ansatz on the form of the exact ground 
state wave function \1> in terms of the unperturbed ground state $ 

l*> = e S \$) , (27) 

where the hermitean operator S is expanded in terms of n-particle and n-hole 
unperturbed states 

S = J2 E ^{k[,...k' n \S n \k 1 ,...k n )a\k' 1 )..J(k' n ) ,a(k n )...a(h)(2S) 

where all the k ' s are hole momenta, i.e. inside the Fermi sphere, and all the k' ' s are 
particle momenta, i.e. outside the Fermi sphere. For translationally invariant systems 
the term S 1 vanishes due to momentum conservation. The exponential form is chosen in 
order to include, as much as possible, only "linked" terms in the expansion of S, in the 
spirit of the linked-cluster theorem discussed above. As already noticed, the unlinked 
diagrams can indeed be summed up by an exponential form. This form also implies 
the normalization ( ( 1 ) | X I / ) = 1. The functions S n are expected to describe the n-body 
correlations in the ground state. As an illustration, let us consider only 62 for simplicity 
and let us assume that it can be considered local in coordinate space, S^r, — r,) = Xij, 
where the labels ij include spin-isospin variables. Then the correlated ground state can 
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be written 

*(ri,r 2 , ) = U t<1 f^( ri ,r 2 , ) , (29) 

where the product runs over all possible distinct pairs of particles and fa = exp(2xij)- 
In general, however, the functions S n are highly non-local in coordinate space and the 
expression for the ground state wave function cannot be written in such a simple form. 

The eigenvalue equation for the exact ground state \1/ can be re-written as a (non- 
hermitean) eigenvalue equation for the unperturbed ground state <£> with a modified 
Hamiltonian, transformed according to a similarity transformation generated by S 

e- § He § \$>) = E\$) . (30) 

The equations for the total energy E and for the correlation functions S n can be 
obtained by multiplying systematically Eq. (30) by the unperturbed ground state, two 
particle-two hole states, three particle-three hole states, and so on. The multiplication 
by ($1 gives a particularly simple expression for the total energy. If only a two-body 
interaction V is present, one gets 

E = (<S>\e~ § H e § \$) = E + (<&\{V + [V, S 2 ]-}\<&) , (31) 

where E is the unperturbed total (kinetic) energy, while all the other terms in the 
expansion of the similarity transformation actually vanish. Therefore, in principle the 
exact total energy can be obtained from the knowledge of the exact two particle- two 
hole amplitude S 2 only. More explicitly 

E = E + X - ]T (hhlWzlhh) (32) 

kl,k2<kp 

where 

(hh^lhh) = (k 1 k 2 \{v + vs 2 }\k 1 k2) 

(33) 

= (fciAfelVlfciAfe) + Eki^kFihk^Vlk'^ik'^lS^k^) . 

Of course, the amplitude S 2 is connected with the higher order amplitudes S 3 ....S n As 

mentioned above, the equations linking the lowest order amplitudes with the higher ones 
are obtained by multiplying Eq. (30) by the unperturbed n particle-n holes bra states (n 
larger or equal to 2). These equations are the constitutive "Coupled Cluster" equations, 
which are equivalent to the eigenvalue equation for the ground state. Approximations 
can be obtained by truncating this chain of equations to a certain order m, i.e. neglecting 
S n for n > m. The meaning of the truncation can be read from the ansatz Eq. (27), it 
amounts to consider correlated n particle n hole components in the ground state up to 
n — m, while higher order components with n > m are just antisymmetrized products 
of the lower ones (note the exponential form, which produces components of arbitrary 
higher orders). 

This form of the CCM equations can be also obtained from the variational principle, 
i.e. by demanding that the mean value of the Hamiltonian in the ground state ^ of Eq. 



Equation of State of Nuclear Matter at high baryon density 



26 



(27) is stationary under an arbitrary variation of the state vector orthogonal to ty. Such 
a variation can be written 

6\V) = e- s1 6Se- s \V) , (34) 

where SS corresponds to an arbitrary variation of the function S n in Eq. (28). It is easily 
verified that such a variation is indeed orthogonal to This is equivalent to take SS 
systematically proportional to a n particle - n hole operator, for all non-zero n, and to 
require that the corresponding energy variation vanishes. This set of conditions gives for 
the functions S n the same CCM equations, which are therefore variational in character. 
The energy can be still taken from Eq. (31), but variants are possible (Navarro 2002). 

However, the CCM equations, as they stand, cannot be applied to calculations in 
nuclear matter or in nuclei. The main correlations in nuclear systems come from the 
strong short range repulsive core, and this part of the NN interaction requires special 
treatment. The many-body wave functions must take into account the overall strong 
repulsion which is present whenever two particles approach each other at a distance 
smaller than the core radius. This requirement must be incorporated systematically 
in the correlation functions S n , otherwise no truncation of the expansion would be 
feasible. The simplest way to proceed is to renormalize the original NN interaction and 
introduce an effective interaction which takes into account the two-body short range 
correlations from the start, so that all the remaining contributions of the expansion are 
expressed in terms of the renormalized, and hopefully "reduced", interaction. In the 
BBG expansion this is done by introducing the G-matrix, and a similar procedure can 
be followed within the CCM scheme. Originally such a line was developed by Kiimmel 
and Luhrmann (1972), and resulted in the so-called Hard Core truncation scheme. A 
similar procedure has been followed recently in finite nuclei (Kowalski 2004 and Dean 
2004), where the G-matrix, first calculated in an extended space, is then used in the 
CCM scheme to calculate systematically the correlations not included in the G-matrix. 
The introduction of the G-matrix of course changes the order of the diagrammatic 
expansion in the CCM method, in particular while the original CCM scheme treats 
particle-particle (short range) and particle-hole correlations (long range) on the same 
footing, the modified CCM scheme shifts the long range part of the correlations to 
higher orders and introduces the G-matrix as the effective two-body interaction in the 
corresponding terms of the expansion. The formal scheme along this lines has been 
developed for nuclear matter by Day (1983). 

Once this procedure is introduced, in the resulting CCM equations the variational 
property is of course lost, as in the case of the BBG expansion. 

Furthermore, a single particle potential U(k) can also be introduced in the CCM 
method. While the original CCM set of equations are formally independent of U (k) at 
each level of truncation, the modified equations do not depend on the single particle 
potential only if no truncation is performed. A detailed analysis of the connection 
between CCM and BBG expansions was presented by Day (1983). In the modified 
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CCM equations, one introduces the effective interaction 

W = - y2(k!k 2 \V I k 3 k 4 )a k la k l (e~ S a k4 a k3 e S ) , (35) 
2 ft} V Jc 

where the operator in parenthesis, if applied to the unperturbed ground state, can 
produce two hole states, 3 holes and 1 particle, 4 holes and two particles and so on. The 
subscript c indicates that, in this expansion, only the terms which do not annihilate the 
unperturbed ground state are retained, i.e. no a\ with k < kp or a k with k > hp are 
retained. The operator W can be also expanded in n particles - n holes operators 

W = J2 E ^(A;' 1 ,...^|^ n |A; 1 ,...A; n )a t (A; , 1 )...at(A;;)a(A; n )...a(A; 1 ) (36) 

in exactly the same fashion as the operator S. The functions W n are related with the 
functions S n . Schematically this relation can be written 

W n = V5 n , 2 + VSn-! + VS n + £ VS k S n - k . (37) 

k<n~2 

The modified CCM equations are obtained by the same procedure as before. The 
equations involves now both the functions W n and the functions S n . Together with the 
relationship Eq. (37), a closed set of equations is then obtained, which is again equivalent 
to the original eigenvalue problem for the ground state. The ground state energy is still 
given by Eq. (32), since the relation between W 2 and S 2 , according to Eq. (37) for n = 2, 
is indeed given by Eq. (33). The truncation scheme (the "Bochum" truncation scheme) 
is now performed on both W n and S n , i.e. the truncation at order m corresponds to 
neglecting the functions W n and S n for n > in. If we truncate the expansion at m — 2, 
only W 2 and S 2 are retained. The quantity W 2 can be readily identified with the on- 
shell G-matrix of the BBG expansion and the function S 2 with the corresponding defect 
function. If the self-consistent single particle potential is introduced, one then gets at 
this level exactly the Brueckner approximation. 

As already discussed, the G-matrix can be introduced in all the terms of the 
Coupled-Cluster expansion. In this case each term of the expansion coincides with 
one diagram in the BBG method. However, it turns out that the ordering of terms 
according to the modified CCM truncation scheme at increasing n does not coincide 
completely with the ordering of diagrams in the hole-line expansion, i.e. for n > 2 the 
CCM expansion at a given truncation n includes also diagrams with a number of hole 
lines larger than n. In particular, the truncation at n = 3 includes also "ring diagrams" 
with an arbitrary number of hole lines, i.e. the whole series of the particle-hole ring 
diagrams initiated by the diagram of Fig. 14, adding more and more particle-hole 
bubbles. These have been shown to be small (Day 1981), and therefore the CCM 
can be considered equivalent to the BBG expansion up to the three hole-line level of 
approximation. 

The Coupled Cluster method gives a new insight into the structure and meaning 
of the hole- line expansion according to the BBG method. In fact, as we have seen the 
CCM is based on the ansatz (27) for the ground state wave function, and it is likely that 
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the same structure of the ground state is underlying the BBG expansion. At Brueckner 
level it is then consistent to assume that the ground state wave function is given by 

|*Br«> = e^ 2 |$) (38) 
with S 2 the Brueckner defect function. 

5. The variational method 

The variational method for the evaluation of the ground state of many-body systems was 
developed since the formulation of quantum theory of atoms and molecules. It acquires 
a particular form in nuclear physics because of the peculiarities of the NN interaction. 
The strong repulsion at short distance has been treated by introducing a Jastrow-like 
trial wave function. The complexity of the NN interaction needs special treatment and 
the introduction of more complex correlation factors. Many excellent review papers exist 
in the literature on the variational method and its extensive use for the determination 
of nuclear matter EoS (Navarro et al. 2002, Pandharipande and Wiringa 1979). Here 
we restrict the exposition to the essential ingredients of the method to the purpose of a 
formal and numerical comparison with the other methods. 

In the simple case of a central interaction the trial ground state wave function is 
written as 

*(ri,r 2 , ) = U^jfinjMn,^, ) , (39) 

where $ is the unperturbed ground state wave function, properly antisymmetrized, and 
the product runs over all possible distinct pairs of particles. The similarity with the 
wave function of Eq. (29) is apparent and indicates a definite link with BBG and CCM 
methods. The correlation function /(r^) is here determined by the variational principle, 
i.e. by imposing that the mean value of the Hamiltonian gets a minimum (or in general 
stationary point) 

5 (y\Hm , . 

— — — — = . (40) 

In principle this is a functional equation for the correlation function /, which however 
can be written explicitly in a closed form only if additional suitable approximations are 
introduced. A practical and much used method is to assume a parametrized form for 
/ and to minimize the energy with respect to the set of parameters which constrain its 
form. Since, as previously discussed, the wave function is expected to decrease strongly 
whenever two particles are at distance smaller than the repulsive core radius of the NN 
interaction, the function /(r^) is assumed to converge to 1 at large distance and to go 
rapidly to zero as — > 0, with a shape similar to the one shown in Fig. 18 for the 
correlated two-body wave function. Furthermore, at distance just above the core radius 
a possible increase of the correlation function beyond the value 1 is possible. 

For nuclear matter it is necessary to introduce a channel dependent correlation 
factor, which is equivalent to assume that / is actually a two-body operator F^. One 
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then assumes that F can be expanded in the same spin-isospin, spin-orbit and tensor 



operators appearing in the NN interaction. Momentum dependent operators, like spin- 
orbit, are usually treated separately. The product in Eq. (39) must be then symmetrized 
since the different terms do not commute anymore. The most flexible assumption 
on the F's is to impose that they go to 1 at a given "healing" distance d with zero 
derivative. The healing distances, which eventually can be defined for each spin-isospin 
and tensor channels, are then taken as variational parameters, while the functions for 
r < d are determined directly from the variational procedure. In principle, the condition 
of energy minimum (or extremal) should produce a set of Euler-Lagrange equations 
which determine the correlation factors. In practice, a viable explicit form can be used 
only for the two-body cluster terms, as discussed below. 

If the two-body NN interaction is local and central, its mean value is directly related 
to the pair distribution function g(r) 



The main job in the variational method is to relate the pair distribution function 
to the correlation factors F. In general this cannot be done exactly, and one has to 
rely on some suitable expansion. For the central part of the correlations, the physical 
quantity which describes the main perturbation with respect to the free Fermi gas is 
the function 1 — F(r) 2 = h(r), which is a measure of the strength of the short range 
part of the correlation. One can then expand the square of the correlated wave function 
in the components with a given number of /i-factors, and correspondingly the energy 
mean value can be expanded in different terms, each one with a given number of h- 
functions. If the full NN interaction is considered, also the non central component 
of the correlation factors, F nc , must be included in the expansion. In this case the 
smallness factors are F% c and the product F nc ■ h, since they are expected to be small 
and vanish at large distance. The different terms can be represented graphically by 
diagrams to help their classification and identify their possible cancellations. It turns 
out (Fantoni and Rosati 1974) that the mean value, at least in the thermodynamic 
limit, is the summation of the so-called "irreducible" diagrams, in strong similarity 
with the linked-cluster theorem of the BBG expansion. Indeed, the reducible diagrams 
are canceled exactly by the expansion of the denominator in Eq. (42). The problem 
of calculating g(r) from the ansatz of Eq. (39) has also a strong similarity with the 
statistical mechanics of a classical gas at finite temperature, where different methods 
to sum up infinite series of diagrams in the so-called "virial expansion" have been 
developed, noticeably the Hypernetted Chain (HNC) summation method (Leeuwven et 
al. 1959). These methods can be almost literally translated to the case of boson systems. 
With some modifications due to the different statistics, they can be extended (Fantoni 
and Rosati 1974) to fermion systems (FHNC), provided in this case the correlations are 




(41) 



where 




(42) 
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taken to be only central ("Jastrow type" correlations) in Eq. (39), i.e. the correlation 
factors are assumed to be only dependent on the coordinates. Unfortunately, in nuclear 
matter, as already mentioned, correlations are of complicated structure due to the NN 
interaction, and the HNC method can be applied only within approximate schemes, 
like the Single Operator Chain (SOC) summation method (Pandharipande and Wiringa 
1979, Lagaris and Pandharipande 1980, Lagaris and Pandharipande 1981), called also 
Variational Summation Method (VSM). In VMS only chains with a given correlation 
operator are considered. In general, the correlation functions are calculated at the two- 
body cluster level, where one gets the Euler-Lagrange coupled equations for all operator 
channels, that can be solved exactly for the set of correlation functions F p (rij) at fixed 
values of the "healing distances" . The index p labels the different two-body operators, 
spin-spin, spin-isospin, tensor, and so on. The VMS method is then applied, keeping 
the same set of correlation functions, to calculate the total energy. The procedure 
is repeated for different values of the healing distances and the energy minimum is 
found within this parameter space. The minimization gives of course automatically 
the ground state wave function. The VSM allows one to include a definite class of 
higher order "clusters" beyond the two-body ones. However, particle clusters in the 
variational method are physically quite different from the ones in the CCM method as 
well as from the BBG ones, where particle "clusters" are defined in terms of diagrams 
with a given number of hole-lines, according to the hole-line expansion. This point 
will be discussed in the next section. Generally speaking, the summation of clusters 
performed by chain summations are expected to include long range correlations, while 
the variational procedure leading to the Euler-Lagrange equations should include mainly 
short range correlations. Indeed, in the low density limit the Euler-Lagrange equation 
reduces to the Schroedinger equation for two particles in free space. 

6. A critical comparison 

As it has been shown by Jackson et al. (1982), in the low density limit, where two- 
body correlations dominate, the Euler-Lagrange equations of the variational method 
are equivalent to the summation of the ladder diagrams of the BBG expansion, while 
the hypernetted chain summation is related to the ring diagram series. This result is 
actually valid only for boson systems, while for Fermi systems it holds approximately, 
only by means of a suitable averaging over entry energy and momenta of the diagrams 
appearing in the BBG expansion. Indeed, the correlation factors are at most state 
dependent in the variational approach, while in principle they should depend also on 
both energy and total momentum. In any case, due to these approximate links, it was 
suggested (Jackson et al. 1982) to use in many-body systems in general a "parquet" 
summation, where both particle-particle short range correlations and chain summations 
of the ring type are treated on the same footing. However, this method has never been 
systematically exploited in the case of nuclear matter, and therefore the approach will 
not be discussed here. 
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The most relevant difference between the BBG (or CCM) method and the 
variational one is the introduction of the self-consistent single particle potential, which 
is not explicitly introduced in the variational procedure. As already noticed, with this 
modification the CCM and the BBG expansion are not any more of variational character, 
in general, at a given level of truncation. However, at the same time a large fraction 
of higher order correlations are effectively embodied in the single particle potential and 
the speed of convergence of the expansion is substantially improved. In the variational 
approach the average single particle potential is implicitly built up along the cluster 
expansion. It is likely that this is the reason of the slow convergence in the order of the 
clusters included in the chain summations (Morales et al. 2002), and also for this reason 
the meaning of "clusters" is not straightforwardly the same in the different methods. 

In the variational approach three-body correlations arise as cyclic products of three 
two-body factors, e.g. f{rij)f{r-jk)f(j'ki)- This contribution has been recently (Morales 
et al. 2002) calculated exactly in symmetric and pure neutron matter for realistic 
interactions. Irreducible three-body correlations can be introduced from the start by 
multiplying the uncorrelated wave function not only by two-body correlation factors 
f(rij) but also by three-body correlation factors fijk, which will then include those 
three-body correlations which cannot be be expressed as product of two-body ones. As 
noticed by Luhrmann (1975), this also indicates a difference with the BBG (and CCM) 
expansion, where the whole three-body correlations are included in the energy term 
generated by the Bethe-Fadeev equations. 

Despite all these differences, some similarity of the methods appear to be present, 
while a more detailed comparison can be made only at the level of the numerical results. 

In summary, the main differences between the variational and the BBG approaches 
can be identified as follows. 

1. In the BBG method for the nuclear EoS the kinetic energy contribution is kept at its 
unperturbed value at all orders of the expansion, while all the correlations are embodied 
in the interaction energy part. This characteristic of the BBG method is not due to any 
approximation but to the expansion method, where the modification of the occupation 
numbers due to correlations is treated on the same footing and at the same order as the 
other correlation effects. In the variational method both kinetic and interaction parts 
are directly modified by the correlation factors. 

2. The correlation factors introduced in the variational method are assumed to be 
essentially local, but usually state dependent. The corresponding (implicit) correlation 
factors in the BBG expansion are in general highly non-local and energy dependent, 
besides being state dependent. 

3. In the BBG method the auxiliary single particle potential U (k) is introduced within 
the expansion in order to improve the rate of convergence. No single particle potential 
is introduced in the variational procedure for the calculation of the ground state energy 
and wave function. Of course, once the variational calculation is performed, the single 
particle potential can be extracted. This also should imply that the rate of convergence 
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in terms of the order of the clusters is slower in the variational method. It was indeed 
shown by Morales et al. (2002) that one needs clusters at least up to order 5 to 
get reasonable convergence, but in principle this does not create any problem in the 
variational method (while in the BBG expansion it would be a disastrous difficulty). It 
has to be stressed anyhow that the physical meaning of "cluster" is quite different in the 
two methods, being more related to long range correlations in the variational scheme, 
to the short range ones in BBG. 

Point 3 is probably the most relevant difference between the two methods, but in 
any case it is difficult to estimate to which extent each one of the listed differences can 
affect the resulting EoS. 

The similarity and connection between the two methods can be found by 
interpreting on physical grounds the diagrammatic expansion used in each one of them. 
The two-body correlations are surely described by the lowest order diagram of the 
variational method, which corresponds to a factor fy, which in turn can be related 
to the G- matrix, i.e. to the Brueckner approximation (with the warning of point 2). 
The hypernetted sums, in their various form, should be connected with the series of 
ring diagrams starting from the one discussed in connection with the three hole-line 
diagrams (including an arbitrary number of loops). As mentioned above, the three-body 
correlations included in the Bethe-Fadeev equations can be related to the irreducible 
product of three fa factors. For boson systems all these connections are more stringent, 
for fermion systems like nuclear matter they are much less transparent and one has to 
rely on physical arguments. 

7. The Equation of State from the BBG and the variational approach 

The first obvious requirement any EoS must satisfy is the reproduction of the so-called 
"saturation point" (SP) , extracted from the fit of the mass formula to the smooth part 
of the binding energy of nuclei along the stability valley. To be definite, we will take 
the values e = —16 MeV and p = 0.17/m™ 3 for the energy per particle and density, 
respectively, as defining the SP of symmetric nuclear matter. As it is well known, no two- 
body force which fits the NN phase shifts was found to be able to reproduce accurately 
the SP. In early applications of Brueckner theory it was realized that the SP predicted 
by different phase-equivalent NN interactions lie inside the so-called "Coester band", 
after Coester et al. (1970). The band misses the phenomenological SP, even taking 
into account the intrinsic uncertainty coming from the extraction procedure (different 
mass formulae, different fit procedures, etc.). The band indicates that either the binding 
energy is too small but the density is correct, or the binding energy is correct but the 
density is too large, see Fig. 23. Furthermore the position along the band was related to 
the strength of the tensor forces, i.e. to the percentage of D-wave in the deuteron. Higher 
values of the strength were corresponding to the upper part of the band. However the 
analysis was done in the gap choice, as discussed above. The use of the continuous choice 
within the Brueckner method changes substantially the results, as shown in Fig. 23. If 
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Figure 23. Saturation curve (full line) for the Argonne V14 potential (Wiringa 
et al. 1984) in the Brueckner approximation with the gap choice for the single 
particle potential. The saturation points for other NN interactions within the same 
approximation arc indicated by open circles. They display the Coester band discussed 
in the text. The reported percentages indicate, for some of the interactions, the 
strength of the D-wave component in the deuteron. Most of the reported systematics 
is reproduced from Machleidt R 1989, Adv. Nucl. Phys. 19, 189, where the force 
corresponding to each label is specified and the corresponding references are given. The 
big square marks the approximate empirical saturation point. The stars correspond to 
some results obtained within the Brueckner scheme with the continuous choice for the 
single particle potential, as discussed in the text. 



some the most modern local forces, with different deuteron D-wave percentages, are used, 
the SP turns out to be restricted in this case to an "island" , which however is still shifted 
with respect to the phenomenological SP. The discrepancy does not appear dramatic. 
Taking into account that the Brueckner approximation is the lowest order in the BBG 
scheme, this result is surely remarkable. Unfortunately, as shown in the previous section, 
higher order contributions, namely the three hole-line diagrams, do not change the 
nuclear EoS appreciably and the discrepancy still persists. The variational method 
gives results in full agreement with this conclusion. The deficiency is evidently not in 
the many-body treatment but in the adopted Hamiltonian. Two possible corrections 
can be devised : many-body forces (to be distinguished from many-body correlations), 
in particular three-body forces, and relativistic effects. As we will mention later, it is 
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well known that the two possible corrections are actually strongly related. Here we will 
consider three-body forces in some detail. 

First we compare in Fig. 24. the BBG and variational EoS (Akmal et al. 1998) 
both for symmetric matter and pure neutron matter without three-body forces in order 
to single out the dependence of the results on the adopted many-body scheme. Since we 
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Figure 24. Symmetric matter (lower curves) and pure neutron matter (upper curves) 
EoS for the Argonne vi 8 NN potential calculated within the BBG (dashed lines) and 
the variational (diamonds) methods. Only two-body forces are included. 



focus on the high density part of the EoS, i.e. above saturation density, the comparison is 
displayed in a wide density range. It has to be stressed that the NN phase shifts constrain 
the NN two-body force up to about 350 MeV in the laboratory, which corresponds to 
a relative momentum of about k\ — 2 fm -1 . Densities corresponding to values of hp 
larger than k\ fall surely in the region where an extrapolation is needed and the NN 
force is untested. For pure neutron matter the agreement between the two theories can 
be considered surprisingly good up to quite high density. For symmetric matter the 
good agreement extends up to about 0.6 fm~ 3 , while at higher density the variational 
EoS is substantially higher than the BBG one. The reason for that is unknown. 

This type of agreement is still present if three-body forces (TBF) are introduced to 
the purpose of getting a SP in agreement with the empirical findings. This can be seen 
in Fig. 25 , where calculations with the Argonne vig interaction and the Urbana model 
for three-body forces are presented. These TBF contain an attractive and repulsive part, 
whose structure is suggested by elementary processes which involve meson exchanges 
and three nucleons but cannot be separated into two distinct nucleon-nucleon interaction 
processes. These TBF are phenomenological in character since the two parameters, 
namely the strength of the attractive and the repulsive terms, cannot be fixed from 
first principles but they are adjusted to reproduce accurately experimental data. In 
Fig. 25 we also report (full line) the EoS of Heiselberg and Hjorth- Jensen (1999), who 
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proposed a modification of the variational EoS of Akmal et al. (1998), which prevents 
its superluminal behaviour at high density (this EoS is in fact softer). 
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Figure 25. Symmetric matter (lower curves) and pure neutron matter (upper curves) 
EoS for the Argonne vi 8 NN potential and three-body forces (TBF), calculated within 
the BBG (dashed lines) and the variational (diamonds) methods. The full lines 
correspond to the modified version of the variational EoS of Hciselbcrg and Hjorth- 
Jensen (1999) 

Few observations are in order. It turns out that the parameters of the three-body 
forces which have been fitted to data on few nucleon systems (triton, 3 He and 4 He) have 
to be modified if the SP has to be reproduced within the phenomenological uncertainty. 
Generally speaking the repulsive part has to be reduced substantially. It could be argued 
that this change of parameters poses a serious problem, since if the TBF model makes 
sense one must keep the same parameters at all densities, otherwise higher order many- 
body forces should be invoked. Actually this is a false problem. In fact the discrepancy 
on the SP cannot be reduced more than few hundred MeV (typically 200-300 KeV), due 
to the intrinsic uncertainty in its position. This discrepancy remains essentially the same 
for few-body systems if one uses the same TBF fitted in nuclear matter, and actually the 
contribution of TBF in few-body systems is quite small. Therefore TBF which allow one 
to describe both few-body systems and the nuclear matter SP do exist, if the accuracy is 
kept at the level of 200-300 KeV in energy per particle, and indeed they are not unique. 
To try a very precise overall fit at the level of 10 KeV or better, as it is now possible in 
the field of few-body systems, appears definitely to be too challenging. There are surely 
higher order terms (four-body forces, retardation effects in the NN interaction, other 
relativistic effects, etc.) which could contribute at this level of precision. In fact the 
contribution of TBF to the energy per particle around saturation is in all cases about 
1-2 MeV, while for few-body systems it is one order of magnitude smaller. Therefore a 
TBF tuned to fit the binding energy of few-body system with an accuracy of 1-10 KeV 
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cannot be extended to fit nuclear matter SP, since this would correspond to a quite 
unbalanced fitting procedure. In other words, 100-200 KeV for the energy per particle 
is the limit of accuracy of the TBF model applied to the nuclear EoS in the considered 
wide range of density. 

Secondly, it has to be noticed that once the SP is reproduced by adjusting the TBF, 
it turns out that the parameters are not the same in the BBG and variational methods, 
i.e. the TBF are not the same. Finally the way of incorporating TBF is simplified in 
the BBG method, namely TBF are reduced to a density dependent two-body force by a 
suitable average over the position and spin-isospin quantum numbers of the third particle 
(Grange et al. 1989). The results presented in Fig. 25 are Brueckner calculations with 
TBF included following this procedure. The agreement between the two curves seems 
to indicate that once the SP is reproduced correctly, the full EoS is determined to a 
large extent up to density as high as 0.6 fm~ 3 . However the conclusion is restricted to 
the particular model for the two-body and three-body forces. The possible dependence 
on the considered forces will be discussed in the next section. 

8. Dependence on the two and three-body forces 

The progress in the accuracy and extension of NN experimental data, as well as in their 
fit by different NN interaction models has been quite impressive. The data range up to 
350 MeV in the laboratory. Going beyond this limit requires the introduction of non- 
nucleonic degrees of freedom (mesons, isobar resonances, etc.) and the corresponding 
inelastic channels. The latter possibility looks quite a complex task, and only few cases 
are present in the literature, notably by Ter Haar and Malfliet (1987). Since in this 
paper we restrict to nucleonic degrees of freedom only, we will consider NN two-body 
interactions which do not include explicitly mesonic or isobaric degrees of freedom. 
Despite that, it is a common paradigm that the NN interaction is determined, at least 
in an effective way, by the exchange of different mesons, and practically all modern NN 
interactions take inspiration for their structure by this assumption, in an explicit or 
implicit way. A really large set of two-body interactions has been developed along the 
years, but nowadays it is mandatory to restrict the possible choices to the most modern 
ones, since they are the only ones which fit the widest and most accurate experimental 
data, on one hand, and are more accurate in the fitting, on the other. One can then 
restrict the set of possible NN two-body interactions to the ones which fit the latest 
data (few thousands of data points) with an accuracy which gives a x 2 /datum close 
to 1. With these requirements, the number of NN interactions reduces quite a bit, 
and only few ones can be considered acceptable. The one which is constructed more 
explicitly from meson exchange processes is the recently developed CD Bonn potential 
by Machleidt (2001), which is the latest one of the Bonn potential series. In principle 
this interaction is the one with the best x 2 value. However the experimental data are 
not always consistent to the needed degree of accuracy and some selection must be done. 
In the same work one can find a detailed analysis of the different data set together with 
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the method of selection which has been followed for the most accurate NN interactions 
and the values of the corresponding x 2 , if one includes the data up to the year 2000. 
Among the most accurate NN interactions one has to include the Argonne vig, already 
discussed. It is constructed by a set of two-body operators which arise naturally in 
meson exchange processes, but the form factors are partly phenomeno logical (except, 
of course, the one-pion exchange). This interaction has been recently modified in the 
1 S'o and 3 Si — 3 D 1 channels with the inclusion of a purely phenomenological short range 
non-local force, which substitutes the original potentials below 1 fm (Doleschall and 
Borbely 2000, Doleschall et al. 2003, Doleschall 2004), usually indicated as IS potential. 
This allows one to reproduce the binding energy of three and four nucleon systems very 
accurately without the inclusion of any TBF, at variance with the original Argonne 
vi§. Also the radii of 3 if and 3 He are accurately reproduced, while the radius of 4 He 
is slightly underestimated (Lazauskas and Carbonell 2004). This potential is phase- 
equivalent to the original interaction, but the off-shell behaviour is modified. Finally 
one can mention the latest potentials of the Nijmegen group (Stocks et al. 1994). They 
have also the ideal value x 2 /datum 1 . However the fit was performed separately 
in each partial wave and the corresponding two-body operator structure cannot have 
the simple form as expected from meson exchange processes. This interaction will be 
discussed only marginally. A larger set of interactions, which includes the old NN 
potentials, can be found in Li K H et al. (2006), where the resulting symmetric nuclear 
matter EoS at BHF level are compared. 

It has to be noticed that the three selected NN interactions, vig, CD Bonn and IS, 
give an increasingly better reproduction of the three-body binding energy and radii, and 
at the same time their non-locality is increasing in the same order. They are all phase 
equivalent to a good accuracy, so that the differences appearing in the nuclear matter 
EoS can be solely due to their different off-shell behaviour. It is well known (Coester 
et al. 1970) that phase equivalent potentials can give a quite different saturation point 
and overall EoS, but here the comparison is restricted to a very definite class of realistic 
accurate NN interactions, with an operatorial structure which is quite similar and 
suggested by meson exchange processes. In Fig. 26 we compare the three corresponding 
EoS (Baldo and Maieron 2005) at the BHF level and with the inclusion of three hole- 
line diagrams (no TBF). First of all one can notice that also for the CD Bonn and IS 
interactions the three hole-line contributions are still relatively small (especially if one 
compares the interaction energies at two and three hole-line levels). The convergence 
of the hole expansion appears to be a general feature. Furthermore, at increasing non- 
locality the EoS become softer and the SP tends to run away from the empirical SP 
and unreasonably large values for the binding energy and the SP density are obtained. 
In the figure two versions of the IS potential are considered, NL1, where only the 1 Sq 
channel is modified, and NL2, where also the 3 Si - 3 Di channel is modified. The results 
indicate that the problem of reproducing the empirical SP with two-body realistic and 
accurate interactions cannot be solved by introducing non-locality, which modifies the 
off-shell properties of the two-body potentials. Furthermore, the requirement of a very 
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Figure 26. Symmetric matter EoS at two hole-line level (upper panel) and at three 
hole-line level (lower panel) for different NN interactions, the Argonnc v 18 , the CD 
Bonn and two versions of the IS potential (NL1 and NL2), see the text for detail. 



The necessity of introducing TBF around saturation is in any case definitely 
confirmed. According to the results shown in Fig. 26. it is apparent that TBF cannot 
be unique, but they depend on the two-body forces employed, since each one of the 
two-body forces gives a different discrepancy for the SP and therefore needs a different 
correction. This is not surprising, since both two-body force and TBF should originate 
within the same physical framework and therefore they are intimately related. In 
particular, in the nucleon-meson coupling models TBF should be generated by processes 
which involve the coupling constants which are already present in the two-body forces. 
Well known examples are depicted in Fig. 27. Other couplings, like meson-meson ones, 
appear only at the TBF level. It has to be stressed that all these processes must be 
considered within an effective theory framework (i.e. a theory with cutoff). The problem 
of consistency between two-body interactions and TBF has been taken systematically 
by Grange et al (1989) and further developed in recent works (Zuo et al. 2002). Since 
processes which include meson-meson couplings seem to be small and can be neglected 
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Fig. 27. Some processes which can contribute to three-nucleon forces. 



in first approximation, TBF calculated along these lines do not contain in principle any 
additional parameters. It is not surprising then that the corresponding EoS has a SP 
which is appreciably more shifted away from the empirical one if compared with the EoS 
with the same two-body force (v 18 ) but with the phenomenological TBF (see Li Z H 
et al. 2006 fe ). In any case TBF can be a good starting point for further improvements. 
It is unclear if other more complex processes can play a role. The effect of these TBF 
in few-body systems is not known, but it is expected to be not very large, since at low 
density also the contribution of these TBF in nuclear matter becomes quite small. The 
most significant difference with the phenomenological TBF is the stiffness of the EoS 
at high density, which turns out to be much higher, as can be seen also in the case of 
pure neutron matter. The energy and pressure rise steeply above saturation, and this 
can create some problems, as discussed later. 

9. The Dirac-Brueckner approach 

As already mentioned, one of the deficiencies of the Hamiltonian considered in the 
previous sections is the use of the non-relativistic limit. The relativistic framework is 
of course the framework where the nuclear EoS should be ultimately based. The best 
relativistic treatment developed so far is the Dirac-Brueckner approach. Excellent review 
papers on the method can be found in the literature (Machleidt 1989) and in textbooks 
(Brockmann and Machleidt 1999). Here we restrict the presentation to the main basic 
elements of the theory and to the latest results, in order to make the comparison with the 
other methods more transparent. We will follow closely the presentation by Brockmann 
and Machleidt (1999) but we will make reference also to the more recent developments. 
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In the relativistic context the only NN potentials which have been developed are 
the ones of OBE (one boson exchange) type. The starting point is the Lagrangian for 
the nucleon-mesons coupling 

C pv = - i^ 7 V^( ps ) (43) 
rrip S 

C s = +9s^ {s) (44) 

A; 
am' 



C v = - g^r^ - ^a^(d^ - dvipf) (45) 



with tp the nucleon and the meson fields, where a indicates the type of meson and 
H the Lorentz component in the case of vector mesons. For isospin 1 mesons, is 
to be replaced by r • f^ a \ with r l (I = 1,2,3) the usual Pauli matrices. The labels 
ps, pv, s, and v denote pseudoscalar, pseudovector, scalar, and vector coupling/field, 
respectively. 

The one-boson-exchange potential (OBEP) is defined as a sum of one-particle- 
exchange amplitudes of certain bosons with given mass and coupling. The main 
difference with respect to the non-relativistic case is the introduction of the Dirac-spinor 
amplitudes. The six non-strange bosons with masses below 1 GeV/c 2 are used. Thus, 

Vobep= E V«° BE (46) 

a=TT,ri,p,LO,S,(T 

with 7r and rj pseudoscalar, a and 5 scalar, and p and uj vector particles. The 
contributions from the isovector bosons 7r, 5 and p contain a factor T\ • r 2 . In the 
so called static limit, i.e. treating the nucleons as infinitely heavy (their energy equals 
the mass) the usual denominator of the interaction amplitude in momentum space, 
coming from the meson propagator, is exactly the same as in the non-relativistic case 
(since in both cases meson kinematics is relativistic). This limit is not taken in the 
relativistic version, noticeably in the series of Bonn potentials, and the full expression 
of the amplitude with the nucleon relativistic (on-shell) energies is included. As an 
example, let us consider one pion exchange. As it is well known, in the non-relativistic 
and static limit the corresponding local potential in momentum space reads (in standard 
notations) 

4M 2 k 2 c 2 + {mc 2 ) 2 1 1 V 

with k = q — q', where q and q' are the initial and final relative momenta of the 
interacting nucleons. This has to be compared with the complete expression of the 
matrix element between nucleonic (positive energy) states (Machleidt 2000). In the 
center of mass frame it reads 

yfuii = 9l (E' + M)(E + M) ( a, ■ q' _ a x ■ q \ ( a 2 ■ q! _ a 2 ■ q \ 
4M 2 k 2 c 2 + (mc 2 ) 2 \E' + M E + MJ \E' + M E + MJ 



where E, E' are the initial and final nucleon energies. One can see that in this case 
some non-locality is present, since the matrix element depends separately on q and q'. 
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Putting E — E' — M, one gets again the local version. Notice that in any case the 
two versions coincide on-shell (E = E'), and therefore the non- locality modifies only 
the off-shell behaviour of the potential. The matrix elements are further implemented 
by form factors at the NN-meson vertices to regularize the potential and to take into 
account the finite size of the nucleons and the mesons. In applications of the DBHF 
method usually one version of the relativistic OBE potential is used, which therefore 
implies that a certain degree of non- locality is present. As already anticipated in the 
previous section, this is also true if these potentials are used within the non-relativistic 
BHF method. 

The fully relativistic analogue of the two-body scattering matrix is the covariant 
Bethe-Salpeter (BS) equation. In place of the NN non-relativistic potential the sum 
V of all connected two-particle irreducible diagrams has to be used, together with the 
relativistic single particle propagators. Explicitly, the BS equation for the covariant 
scattering matrix T in an arbitrary frame can be written 



7V, q\P) = V(q', q\P) + f d 4 kV(q', k\P)G{k\P)T{k, q\P) , (47) 



with 



G(k\P) 



1 



1 



(2tt) 4 a r- 



(2k) 



fl-M + ie)W (± P 



[\P + k) 2 — M 2 + ie 



2 

(i) 



-M + ie)( 2 ) 
\ P- 16 + M 



[\P — k) 2 — M 2 + ie 



(2) 



(48) 



(49) 



where q, k, and q' are the initial, intermediate, and final relative four-momenta, 
respectively (with e. g. k — (/c ,k)), and P = (P ,P) is the total four-momentum; 
jk = 7 M A; M . The superscripts refer to particle (1) and (2). Of course all quantities 
are appropriate matrices in spin (or helicity) and isospin indices. The use of the OBE 
potential as the kernel V is equivalent to the so-called ladder approximation, where one 
meson exchanges occur in disjoint time intervals with respect to each other, i.e. at any 
time only one meson is present. Unfortunately, even in the ladder approximation the 
BS equation is difficult to solve since V is in general non-local in time, or equivalently 
energy dependent, which means that the integral equation is four-dimensional. It is 
even not sure in general if it admits solutions. It is then customary to reduce the four- 
dimensional integral equation to a three-dimensional one by approximating properly the 
energy dependence of the kernel. In most methods the energy exchange k is fixed to 
zero and the resulting reduced BS equation is similar to its non-relativistic counterpart. 
In the Thompson reduction scheme this equation for matrix elements between positive- 
energy spinors (cm. frame) reads 

r d 3 k M 2 1 

T«, q) = V iq! , q) + / ^VM, k) ^ ^-^_T(k, q|P) (50) 

where both V(q', q) and T have to be considered as matrices acting on the two-particle 
helicity (or spin) space, and _E k = Vk 2 + M 2 is the relativistic particle energy. In 
the alternative Blankenbecler-Sugar (Machleidt 2000) reduction scheme some different 
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relativistic kinematical factors appear in the kernel. This shows that the reduction 
is not unique. The partial wave expansion of the T-matrix can then be performed 
starting from the helicity representation. The corresponding amplitudes include single 
as well as coupled channels, with the same classification in quantum numbers JLS as 
in the non relativistic case and therefore their connection with phase shifts is the same 
(Brockmann and Machleidt 1998). In the intermediate states of momentum k only the 
positive energy states are usually considered (by the proper Dirac projection operator). 
As in the case of the OBEP potential, again the main difference with respect to the 
non-relativistic case is the use of the Dirac spinors. 

The DBHF method can be developed in analogy with the non-relativistic case. The 
two-body correlations are described by introducing the in-medium relativistic G-matrix. 
The DBHF scheme can be formulated as a self-consistent problem between the single 
particle self-energy £ and the G-matrix. Schematically, the equations can be written 

G = V + i J VQggG 

£ = -i JjTr[gG] — gG) (51) 

where Q is the Pauli operator which projects the intermediate two particle momenta 
outside the Fermi sphere, as in the BHF G-matrix equation, and g is the single particle 
Green' s function. The self consistency is entailed by the Dyson equation 

9 = 9o + 9o^9 

where go is the (relativistic) single particle Green's function for a free gas of nucleons. 
The self-energy is a matrix in spinor indices, and therefore in general it can be expanded 
in the covariant form 

E(fc, kp) = E s (k, kp) - 7o S (A;, k F ) + 1 ■ kE„ (52) 

where 7 M are the Dirac gamma matrices and the coefficients of the expansion are scalar 
functions, which in general depend on the modulus |k| of the three-momentum and on 
the energy fc - Of course they also depend on the density, i.e. on the Fermi momentum 
k F - The free single particle eigenstates, which determine the spectral representation of 
the free Green' s function, are solutions of the Dirac equation 

[ 7„Jfc" - M } u(k) = 

where u is the Dirac spinor at four-momentum k. For the full single particle Green's 
function g the corresponding eigenstates satisfy 

[ 7^ - M + E ] u(k)* = 

Inserting the above general expression for £, after a little manipulation, one gets 



[ 7^* - M* }u(k)* = 
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with 



k [ 



,o* 



fc° + S, 
1 + Sf 







fc** = fc* 



M* = 



M + E, 



s 



(53) 



This is the Dirac equation for a single particle in the medium, and the corresponding 
solution is the spinor 



In line with the Brueckner scheme, within the BBG expansion, in the self-energy of 
Eq. (51) only the contribution of the single particle Green' s function pole is considered 
(with strength equal one). Furthermore, negative energy states are neglected and one 
gets the usual self-consistent condition between self-energy and scattering G-matrix. 
The functions to be determined are in this case the three scalar functions appearing in 
Eq. (52). However, to simplify the calculations these functions are often replaced by 
their value at the Fermi momentum. 

In any case, the medium effect on the spinor of Eq. (54) is to replace the vacuum 
value of the nucleon mass and three-momentum with the in-medium values of Eq. 
(53). This means that the in-medium Dirac spinor is "rotated" with respect to the 
corresponding one in vacuum, and a positive (particle) energy state in the medium has 
some non-zero component on the negative (anti-particle) energy state in vacuum. In 
terms of vacuum single nucleon states, the nuclear medium produces automatically anti- 
nucleon states which contribute to the self-energy and to the total energy of the system. 
It has been shown by Brown et al. (1987) that this relativistic effect is equivalent to 
the introduction of well defined TBF at the non-relativistic level. These TBF turn out 
to be repulsive and consequently produce a saturating effect. The DBHF gives indeed 
in general a better SP than BHF. Of course one can wonder why these particular TBF 
should be selected, but anyhow a definite link between DBHF and BHF + TBF is, in 
this way, established. Indeed, including in BHF only these particular TBF one gets 
results close to DBHF calculations, see e.g. Li K H (2006). 

Despite the DBHF is similar to the non-relativistic BHF, some features of this 
method are still controversial. The results depend strongly on the method used to 
determine the covariant structure of the in-medium G-matrix, which is not unique 
since only the positive energy states must be included. It has to be stressed that, in 
general, the self-energy is better calculated in the matter reference frame, while the 
G-matrix is more naturally calculated in the center of mass of the two interacting 
nucleons. This implies that the G-matrix has to be Lorentz transformed from one 
reference frame to the other, and its covariant structure is then crucial. Formally, the 
most accurate method appears to be the subtraction scheme of Gross-Boelting et al. 
(1999). Generally speaking, the EoS calculated within the DBHF method turn out to 
be stiffer above saturation than the ones calculated from the BHF + TBF method. 




(54) 
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10. The compressibility at and above saturation density 

Despite several uncertainties, the compressibility of nuclear matter at saturation can be 
considered as determined within a relatively small range, approximately between 220 and 
270 MeV. The different relativistic and non-relativistic microscopic EoS, considered in 
the previous section, turn out to be compatible with these values of the compressibility, 
As we have seen, however, the compressibility (i.e. stiffness) can differ at high enough 
density, which can be relevant for many phenomenological data. 

In heavy ion collisions at intermediate energy, nuclear matter is expected to be 
compressed and to reach densities few times larger than the saturation value. Several 
observable quantities have been devised that should be sensitive to the stiffness of 
the nuclear EoS. In particular the measure of different types of "flow" is considered 
particularly useful and this line has been followed in many experiments. A more 
ambitious, and probably more questionable, analysis was performed by Daneliewicz 
et al. (2002). The authors consider both the in-plane transverse flow and the elliptic 
flow measured in different experiment on Au + Au collision at energies between 0.2 
and fO GeV/A. According to relativistic Boltzmann-Uehling-Uhlenbeck simulations it 
is claimed that density up to 7 times saturation density is reached during the collisions 
(at the highest energy), and from the data an estimate of the pressure is extracted. 
Together with an evaluation of the uncertainty, the analysis results in the determination 
of a region in the pressure-density plane where the nuclear EoS should be located. In 
this way it appears easy to test a given EoS, and it became popular to confront the 
various microscopic and phenomenological EoS with this region, which is assumed to be 
the allowed one (essentially for symmetric nuclear matter). If one believes the validity 
of this analysis, it turns out that the test is quite stringent, despite the fact that in the 
same work it is also shown that the value of the compressibility at saturation is not at 
all well determined. This means that also at the phenomenological level the value of 
the compressibility at saturation does not determine the EoS stiffness at high density. 
In Fig. 28 the set of microscopic EoS already discussed are reported in comparison with 
the allowed region. The variational EoS of Akmal et al. (1998), as well as the EoS 
derived from BHF together with phenomenological TBF, look in agreement with the 
phenomenological analysis, while the EoS from the DBHF calculations of van Dalen et 
al. (2005) is only marginally compatible, see also the paper by Klahn et al. (2006). 
The non-relativistic EoS calculated with BHF and "ab initio" TBF reported by Zuo et 
al. (2002) looks close to the BHF results with phenomenological TBF up to 2-3 time 
saturation density, giving further support to the phenomenological TBF, see also Zhou 
et al. (2006). However, at higher density it becomes too stiff and definitely falls outside 
the allowed region. 

In turns out that also many phenomenological EoS do not pass this test (Klahn et 
al. 2006). 

It has to be noticed that the flow values, and other observable quantities, in 
general do not depend only on the nuclear EoS, as embodied in the single particle 
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Figure 28. Different EoS in comparison with the phcnomcnological constraint 
extracted by Danielewicz et al. (2002) (shaded area), where p = 0.16 fm -3 . Full line: 
EoS from the BBG method with phenomenological TBF (Zuo et al. 2004). Dashed 
line : modified variational EoS of Heiselberg and Hjorth- Jensen (1999). Dotted line : 
variational EoS of Akmal et al. (1998). Open circle : EoS from the BBG method with 
"ab initio" TBF (Li 2006°). Dash-dotted line : EoS from Dirac-Brueckncr method 
(van Dalcn et al. 2005). 



potential, but also on the in-medium nucleon-nucleon collision cross section and on 
the effective mass (i.e. momentum dependence of the single particle potential). The 
extraction of meaningful information from the experiments requires a careful analysis 
and interpretation of the data. Other quantities which are related to the EoS are the 
rates of particle production, in particular K + and K~ and their ratio. In fact the strange 
particle production is intimately related to the density reached during the collision and 
therefore to the nuclear matter compressibility. However in order to reach meaningful 
conclusions, it is necessary to have a reasonably good description of the behaviour of 
kaons in the nuclear medium at high density, which is not an easy task theoretically. 

On the astrophysical side, as it is well known, each EoS for asymmetric matter 
gives rise to a definite relationship between the mass and radius of neutron stars (NS). 
This is because ordinary NS are bound by gravity and the solution of the Tolmann- 
Opennheimer-Volkoff (TOV) equation, based only on general relativity and the adopted 
EoS, provides the full density profile of the star. Unfortunately it is quite difficult to 
get from observation both the mass and the radius of a single NS, and up to now 
the accuracy, especially of the radius value, is not good enough to discriminate among 
different EoS. The quantity which has been mostly given attention is then the maximum 
mass of NS. The mass vs. radius plot has indeed a maximum value at the smaller radius, 
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beyond which the star configuration is unstable towards collapse to a black hole. This 
maximum is characteristic of each EoS, and the observation of a NS with a mass larger 
than the predicted one would rule out the corresponding EoS. If one assumes that only 
nucleonic degrees of freedom are present inside the NS, then one can adopt the EoS 
discussed above. It turns out that BHF EoS give a maximum mass close to two solar 
masses, while the DBHF and variational ones a slightly larger value, around 2.2 - 2.3 
solar masses. However, as already mentioned, the variational EoS of Akmal et al. (1998) 
becomes superluminal at high density, and actually the corrected one by Heiselberg and 
Hjorth- Jensen (1999) gives a maximum mass very close to 2 solar masses. Also the BHF 
EoS with the TBF by Zuo et al. (2002) gives a maximum mass larger than 2 (Zhou et 
al. 2004), but unfortunately this EoS becomes also superluminal already at relatively 
low density. From the astrophysical observations (mainly binary systems) the masses 
of NS were found, up to few years ago, mostly concentrated around 1.5 solar masses, 
the most precise one being 1.44 (Hulse and Taylor 1975). These values look compatible 
with the theoretical predictions. However, the situation for the maximum mass is by 
far much more complex. First of all it is likely that other degrees of freedom, besides 
the nucleonic ones, can appear inside a NS, in particular hyperonic matter. The BBG 
scheme has been extended to matter containing hyperons in several papers. If the most 
recent hyperon-nucleon and hyperon-hyperon interactions are used the maximum mass 
of NS drops to values below the observational limit (Schulze et al. 2006). There are two 
possibilities to overcome this failure in reproducing the observational constraint. The 
hyperon-nucleon and hyperon-hyperon interactions are poorly known from laboratory 
experiments, which presently are able to provide only few data points to be fitted. 
Different interactions, still compatible with phenomenology, could be able to produce 
a stiffer EoS at high density and consequently larger mass values. Another possibility 
is the appearance of other degrees of freedom, in particular the transition to quark 
matter could occur in the core of NS. This possibility has been studied extensively by 
many authors and indeed it has been found that the onset of the deconfined phase is 
able to increase the maximum mass to values compatible with the observational limit 
and ranging from 1.5 to 1.8 solar masses. These results have been obtained within 
simple quark matter models, like the MIT bag model, the Color Dielectric Model and 
the Nambu-Jona Lasinio model, with the possible inclusion of color superconductivity 
(Drago et al. 2005). If perturbative-like corrections to the simple MIT model are 
introduced, masses up to about 2 can be obtained (Alford et al. 2005). In any case, all 
that shows clearly the great relevance of the nowadays standard observations on NS to 
our knowledge of the high density nuclear EoS. The astrophysical observations are able 
to rule out definite EoS or put constraints on them. 

As a final remark on this subject one has to mention the recent claims of the 
observation of NS with mass definitely larger than 2 (Nice et al. 2005, Ozel 2006). 
If confirmed, these observations would put serious constraints on the nuclear EoS and 
would point to an additional repulsion which should be present at high density, i.e. a 
larger stiffness of the quark matter EoS. It has to be stressed that the nuclear EoS 
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appropriate to NS cannot be directly applied to heavy ion collisions. The NS matter 
is in beta equilibrium and the strange content is determined by chemical equilibrium, 
which cannot be established during the collision time of heavy ions. In fact the hyperon 
multiplicity in heavy ion collisions is much smaller than one and no strange matter can 
be actually formed. Furthermore, the asymmetry of NS matter is much larger than the 
values reachable in laboratory experiments. Of course a good microscopic theory must 
be able to connect the two different physical situations within the same many-body 
scheme, which is one of the main challenges of nuclear physics. 

11. Symmetry energy above saturation 

At sub-saturation density the symmetry energy of nuclear matter seems to be under 
control from the theoretical point of view, since the different microscopic calculations 
agree among each other and the results look only marginally dependent on the adopted 
nuclear interaction. The symmetry energies calculated within the BBG scheme (Baldo 
et al. 2004), for different NN interactions (TBF have a negligible effect) are in 
agreement with each other and are reasonably well reproduced by some of the most 
used phenomenological Skyrme forces. Variational or DBHF calculations give very 
similar results. The approximate agreement of phenomenological calculations with the 
microscopic ones does not hold for all Skyrme forces, as shown e.g. by Chen et al. 
(2006), and a wide spread of values is actually found below saturation. The microscopic 
symmetry energy C sym below saturation density can be approximately described by 

C sym = 31.3(p/po) a6 , 

where po is the saturation density. Notice that the exponent is close to the one for a free 
Fermi gas (of course the absolute values are quite different by approximately a factor 
2). 

The symmetry energy at density above saturation can be studied with heavy ion 
reactions in central or semi-central collisions where nuclear matter can be compressed. 
Particle emissions and productions are among the processes which have been widely 
used to this purpose. Generally speaking the signal coming from these studies are 
weak because several competing effects are very often present at the same time and 
they largely cancel out among each other. In the paper by Ferini et al. (2006) the ratio 
between K + and K° rates and between n~ and ir + has been studied through simulations 
of Au + Au central collisions in the energy range 0.8 - 1.8 GeV. These ratios seem to 
be dependent on the strength of the isovector part of the single particle potentials, but 
the dependence is not so strong, due to the compensation between symmetry potential 
effects and threshold effects. In any case it has to be stressed that the behaviour of K 
mesons, or even n mesons, in nuclear matter is a complex many-body problem, which 
complicates the interpretation of the experimental data. 

To this respect one has to notice that it was suggested by Li et al. (1997) that in 
NS the onset of a kaon condensate could be possible. This can happen due to the steep 
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increase inside the star of the electron chemical potential which can finally equal the 
in-medium mass of K~ mass. Since this condensate produces a substantial softening 
of the EoS, the NS maximum mass turns out to be limited to about 1.5 solar masses. 
The possibility of kaon condensation was re-examined recently by Li et al (2006 a ) . In 
any case, this value looks in contradiction with the latest observational data and shows 
once more the great value of the astrophysical studies on NS for our knowledge of dense 
nuclear matter. 

Another process in NS which is sensitive to symmetry energy is cooling. The main 
mechanism of cooling is the direct Urea (DU) process 

n — > p + e~ + 17 e ; p + e~ — > n + v e 

where neutrinos and antineutrinos escape from the star, cooling the object with a time 
scale of the order of million years. Since the chemical potentials of neutrons and protons 
are quite different because of the large asymmetry of the NS matter, the conservation 
of energy and momenta forbid these reactions when the percentage of protons is below 
14% (when muons are also included). The percentage of protons is directly determined 
by the symmetry energy, and therefore the density at which the threshold for DU occurs 
is directly determined by the density dependence of the symmetry energy. At density 
above saturation this threshold can be different for different EoS. In some case, as for 
the EoS of Akmal et al. (1998), it is practically absent up to almost the central density 
of NS, even for the largest masses. In this case other processes, like the indirect Urea 
process, are the dominant cooling mechanism, which are however much less efficient than 
the DU process. Models of NS which do not include the DU process are only marginally 
successful in reproducing cooling data on NS (Yakovlev and Pethick 2004). If the DU 
threshold is at too "low" density, the cooling process can be too fast, even with the 
inclusion of nuclear matter pairing (Klahn 2006), which hinders the DU process. This is 
what occurs for the EoS derived from the DBHF method, whose symmetry energy rises 
steeply with density. The EoS from BHF, with the inclusion of phenomenological TBF, 
give a threshold density for DU process intermediate between these two cases and seem 
to be compatible with cooling data. For EoS with a similar behaviour of the symmetry 
energy the scenario of NS cooling involves a slow cooling for low masses and a fast one 
for the higher masses. Of course a more detailed description of NS cooling requires the 
many-body treatment of the different processes which can contribute (Blaschke et al. 
2004). Finally, the possible onset of quark matter could again change the whole cooling 
scenario, which has then to be reconsidered, but the above general considerations can 
be still applied. 

12. Conclusions 

In this topical review we have presented the microscopic many-body theories, developed 
along the years, on the nuclear Equation of State, where only nucleonic degrees 
of freedom are considered. The results of different approaches have been critically 
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compared, both at the formal and the numerical levels, with special emphasis on 
the high baryon density regime. The non-relativistic Bethe-Brueckner-Goldstone and 
variational methods are in fair agreement up to 5-6 times saturation density. When 
three-body forces are introduced, as required by phenomenology, some discrepancy 
appears for symmetric nuclear matter above 0.6 fm~ 3 . The dependence of the results 
on the adopted realistic two-body forces and on the choice of the three-body forces 
has been analyzed in detail. It is found that a very precise reproduction of data on 
three- and four-body nuclear systems, as well as of the nuclear matter saturation point, 
is too demanding for the present day nuclear force models. In particular, if a very 
accurate description of few body systems is achieved by a suitable off-shell adjustment 
of the two-body forces, i.e. by introducing a non-local component, then the saturation 
point cannot be reproduced with a reasonable precision. However, local forces with 
phenomenological three-body forces are able to give an approximate reproduction both 
of the properties of few-body systems and of the nuclear matter saturation point, with 
a discrepancy on the binding energy per particle below 200-300 KeV. 

The relativistic Dirac-Brueckner approach, as applied to nuclear matter, gives an 
Equation of State above saturation which is stiffer than the non-relativistic approaches. 
Some ambiguities related to the three-dimensional reduction of the fully relativistic 
two-body scattering matrix in the medium have still to be resolved. 

We have then briefly reviewed the observational data on neutron stars and the 
experimental results on heavy ion collisions at intermediate and relativistic energies 
that could constrain the nuclear Equation of State at high baryon density. The 
EoS of the relativistic Dirac-Brueckner approach seems to present some discrepancies 
in comparison with constraints coming form heavy ion collisions and neutron stars 
data. The microscopic non-relativistic EoS turn out to be compatible with the 
phenomenological constraints available up to now. It looks likely that future 
developments in astrophysical observations and in laboratory experiments on heavy ion 
collisions will further constrain the nuclear EoS and give further hints on our knowledge 
of the fundamental processes which determine the behaviour of nuclear matter at high 
baryon density. 
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